PLaSK library
Loading...
Searching...
No Matches
expansion2d.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>
15#include <boost/range/adaptor/reversed.hpp>
16#include <boost/range/adaptor/transformed.hpp>
17using boost::algorithm::clamp;
18
19#include "expansion2d.hpp"
20#include "solver2d.hpp"
21#include "../meshadapter.hpp"
22
23#define SOLVER static_cast<FourierSolver2D*>(solver)
24
25namespace plask { namespace optical { namespace modal {
26
27ExpansionPW2D::ExpansionPW2D(FourierSolver2D* solver): Expansion(solver), initialized(false),
28 symmetry(E_UNSPECIFIED), polarization(E_UNSPECIFIED) {}
29
30void ExpansionPW2D::setPolarization(Component pol) {
31 if (pol != polarization) {
32 solver->clearFields();
33 if (!periodic && polarization == E_TRAN) {
34 polarization = pol;
35 if (initialized) {
36 reset();
37 init();
38 }
39 solver->recompute_integrals = true;
40 } else if (polarization != E_UNSPECIFIED) {
41 polarization = pol;
42 solver->recompute_integrals = true;
43 } else {
44 polarization = pol;
45 }
46 }
47}
48
49
50void ExpansionPW2D::init()
51{
52 auto geometry = SOLVER->getGeometry();
53
54 shared_ptr<MeshAxis> xmesh;
55
56 periodic = geometry->isPeriodic(Geometry2DCartesian::DIRECTION_TRAN);
57
58 auto bbox = geometry->getChild()->getBoundingBox();
59 left = bbox.lower[0];
60 right = bbox.upper[0];
61
62 size_t refine = SOLVER->refine, nrN;
63 if (refine == 0) refine = 1;
64
65 if (symmetry != E_UNSPECIFIED && !geometry->isSymmetric(Geometry2DCartesian::DIRECTION_TRAN))
66 throw BadInput(solver->getId(), "symmetry not allowed for asymmetric structure");
67
68 if (geometry->isSymmetric(Geometry2DCartesian::DIRECTION_TRAN)) {
69 if (right <= 0.) {
70 left = -left; right = -right;
71 std::swap(left, right);
72 }
73 if (left != 0.) throw BadMesh(SOLVER->getId(), "symmetric geometry must have one of its sides at symmetry axis");
74 if (!symmetric()) left = -right;
75 }
76
77 if (SOLVER->ftt == FourierSolver2D::FOURIER_ANALYTIC) {
78 if (!SOLVER->mesh) SOLVER->setSimpleMesh();
79 if (SOLVER->mesh->size() < 2) throw BadInput(SOLVER->getId(), "mesh needs at least two points");
80 if (!geometry->isSymmetric(Geometry2DCartesian::DIRECTION_TRAN) || SOLVER->mesh->at(0) < 0.) {
81 original_mesh = SOLVER->mesh;
82 } else {
83 shared_ptr<OrderedAxis> new_mesh = make_shared<OrderedAxis>(*SOLVER->mesh);
84 original_mesh = new_mesh;
85 auto negate = [](double x) { return -x; };
86 auto transformed = (*original_mesh) | boost::adaptors::reversed | boost::adaptors::transformed(negate);
87 new_mesh->addOrderedPoints(transformed.begin(), transformed.end(), new_mesh->size());
88 }
89 if (!is_zero(original_mesh->at(0) - (symmetric()? -right : left)))
90 throw BadInput(SOLVER->getId(), "first mesh point ({}) must match left geometry boundary ({})",
91 original_mesh->at(0), symmetric()? -right : left);
92 if (!is_zero(original_mesh->at(original_mesh->size()-1) - right))
93 throw BadInput(SOLVER->getId(), "last mesh point ({}) must match right geometry boundary ({})",
94 original_mesh->at(original_mesh->size()-1), right);
95 }
96
97 if (!periodic) {
98 // Add PMLs
99 if (!symmetric()) left -= SOLVER->pml.size + SOLVER->pml.dist;
100 right += SOLVER->pml.size + SOLVER->pml.dist;
101 }
102
103 double L;
104 // N = 3 nN = 5 refine = 5 nrN = 25
105 if (!symmetric()) { // . . 0 . . . . 1 . . . . 2 . . . . 3 . . . . 4 . .
106 L = right - left; // ^ ^ ^ ^ ^
107 N = 2 * SOLVER->getSize() + 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|
108 nN = 4 * SOLVER->getSize() + 1; // N = 3 nN = 5 refine = 4 nrN = 20
109 nrN = refine * nN; // . . 0 . . . 1 . . . 2 . . . 3 . . . 4 . . . 0
110 double dx = 0.5 * L * double(refine-1) / double(nrN); // ^ ^ ^ ^
111 if (SOLVER->ftt == FourierSolver2D::FOURIER_DISCRETE) // |0 1 2 3|4 5 6 7|8 9 0 1|2 3 4 5|6 7 8 9|
112 xmesh = plask::make_shared<RegularAxis>(left-dx, right-dx-L/double(nrN), nrN);
113 } else {
114 L = 2. * right;
115 N = SOLVER->getSize() + 1;
116 nN = 2 * SOLVER->getSize() + 1;
117 nrN = refine * nN; // N = 3 nN = 5 refine = 4 nrN = 20
118 if (SOLVER->dct2()) { // # . 0 . # . 1 . # . 2 . # . 3 . # . 4 . # . 4 .
119 double dx = 0.25 * L / double(nrN); // ^ ^ ^ ^
120 if (SOLVER->ftt == FourierSolver2D::FOURIER_DISCRETE) // |0 1 2 3|4 5 6 7|8 9 0 1|2 3 4 5|6 7 8 9|
121 xmesh = plask::make_shared<RegularAxis>(dx, right-dx, nrN);
122 } else {
123 size_t nNa = 4 * SOLVER->getSize() + 1;
124 double dx = 0.5 * L * double(refine-1) / double(refine*nNa);
125 if (SOLVER->ftt == FourierSolver2D::FOURIER_DISCRETE)
126 xmesh = plask::make_shared<RegularAxis>(-dx, right+dx, nrN);
127 }
128 }
129 size_t M = matrixSize();
130
131 SOLVER->writelog(LOG_DETAIL, "Creating{2}{3} expansion with {0} plane-waves (matrix size: {1})",
132 N, M, symmetric()?" symmetric":"", separated()?" separated":"");
133
134 if (symmetric()) SOLVER->writelog(LOG_DETAIL, "Symmetry is {0}", (symmetry== E_TRAN)? "Etran" : "Elong");
135
136 if (SOLVER->ftt == FourierSolver2D::FOURIER_DISCRETE) {
137 matFFT = FFT::Forward1D(1, int(nN),
138 symmetric()?
139 SOLVER->dct2()? FFT::SYMMETRY_EVEN_2 : FFT::SYMMETRY_EVEN_1 :
140 FFT::SYMMETRY_NONE);
141 } else {
142 if (!periodic) {
143 xmesh = plask::make_shared<OrderedAxis>(*original_mesh->getMidpointAxis());
144 static_pointer_cast<OrderedAxis>(xmesh)->addPoint(SOLVER->mesh->at(0) - 2.*OrderedAxis::MIN_DISTANCE);
145 static_pointer_cast<OrderedAxis>(xmesh)->addPoint(SOLVER->mesh->at(SOLVER->mesh->size()-1) + 2.*OrderedAxis::MIN_DISTANCE);
146 } else
147 xmesh = original_mesh->getMidpointAxis();
148 }
149
150 // Compute permeability coefficients
151 if (periodic) {
152 mag.reset(nN, 0.);
153 mag[0] = 1.;
154 if (polarization != E_TRAN) {
155 rmag.reset(nN, 0.);
156 rmag[0] = 1.;
157 }
158 } else {
159 if (SOLVER->ftt == FourierSolver2D::FOURIER_DISCRETE) {
160 mag.reset(nN, 0.);
161 if (polarization != E_TRAN) rmag.reset(nN, 0.);
162 // Add PMLs
163 SOLVER->writelog(LOG_DETAIL, "Adding side PMLs (total structure width: {0}um)", L);
164 double pl = left + SOLVER->pml.size, pr = right - SOLVER->pml.size;
165 if (symmetric()) pil = 0;
166 else pil = std::lower_bound(xmesh->begin(), xmesh->end(), pl) - xmesh->begin();
167 pir = std::lower_bound(xmesh->begin(), xmesh->end(), pr) - xmesh->begin();
168 for (size_t i = 0; i != nN; ++i) {
169 for (size_t j = refine*i, end = refine*(i+1); j != end; ++j) {
170 dcomplex sy = 1.;
171 if (j < pil) {
172 double h = (pl - xmesh->at(j)) / SOLVER->pml.size;
173 sy = 1. + (SOLVER->pml.factor-1.)*pow(h, SOLVER->pml.order);
174 } else if (j > pir) {
175 double h = (xmesh->at(j) - pr) / SOLVER->pml.size;
176 sy = 1. + (SOLVER->pml.factor-1.)*pow(h, SOLVER->pml.order);
177 }
178 mag[i] += sy;
179 if (polarization != E_TRAN) rmag[i] += 1. / sy;
180 }
181 mag[i] /= (double)refine;
182 if (polarization != E_TRAN) rmag[i] /= (double)refine;
183 }
184 // Compute FFT
185 FFT::Forward1D fft(1, int(nN),
186 symmetric()?
187 SOLVER->dct2()? FFT::SYMMETRY_EVEN_2 : FFT::SYMMETRY_EVEN_1 :
188 FFT::SYMMETRY_NONE);
189 fft.execute(mag.data());
190 if (polarization != E_TRAN) fft.execute(rmag.data());
191 } else {
192 throw NotImplemented(SOLVER->getId(), "analytic Fourier transform for non-periodic structure"); //TODO
193 }
194 // Smooth coefficients
195 if (SOLVER->smooth) {
196 double bb4 = PI / L; bb4 *= bb4; // (2π/L)² / 4
197 for (std::size_t i = 0; i != nN; ++i) {
198 int k = int(i); if (!symmetric() && k > int(nN/2)) k -= int(nN);
199 mag[i] *= exp(-SOLVER->smooth * bb4 * k * k);
200 if (polarization != E_TRAN) rmag[i] *= exp(-SOLVER->smooth * bb4 * k * k);
201 }
202 }
203 TempMatrix temp = getTempMatrix();
204 cmatrix work(temp);
205 cmatrix workij(N, N, work.data());
206 make_permeability_matrices(workij);
207 }
208
209 // Allocate memory for expansion coefficients
210 size_t nlayers = solver->lcount;
211 coeffs.resize(nlayers);
212 diagonals.assign(nlayers, false);
213 coeff_matrices.resize(nlayers);
214
215 mesh = plask::make_shared<RectangularMesh<2>>(xmesh, solver->verts, RectangularMesh<2>::ORDER_01);
216
217 initialized = true;
218}
219
220void ExpansionPW2D::reset() {
221 coeffs.clear();
222 coeff_matrices.clear();
223 coeff_matrix_mxx.reset();
224 coeff_matrix_rmyy.reset();
225 initialized = false;
226 mesh.reset();
227 mag.reset();
228 rmag.reset();
229 temporary.reset();
230}
231
232
233void ExpansionPW2D::beforeLayersIntegrals(dcomplex lam, dcomplex glam) {
234 SOLVER->prepareExpansionIntegrals(this, mesh, lam, glam);
235}
236
237
238void ExpansionPW2D::layerIntegrals(size_t layer, double lam, double glam)
239{
240 auto geometry = SOLVER->getGeometry();
241
242 double maty;
243 for (size_t i = 0; i != solver->stack.size(); ++i) {
244 if (solver->stack[i] == layer) {
245 maty = solver->verts->at(i);
246 break;
247 }
248 }
249
250 double L = (right-left) * (symmetric()? 2. : 1.);
251
252 bool epsilon_isotropic = true, epsilon_diagonal = true;
253
254 if (SOLVER->ftt == FourierSolver2D::FOURIER_DISCRETE) {
255 size_t refine = SOLVER->refine;
256 if (refine == 0) refine = 1;
257
258 #if defined(OPENMP_FOUND) // && !defined(NDEBUG)
259 SOLVER->writelog(LOG_DETAIL, "Getting refractive indices for layer {}/{} (sampled at {} points) in thread {}",
260 layer, solver->lcount, refine * nN, omp_get_thread_num());
261 #else
262 SOLVER->writelog(LOG_DETAIL, "Getting refractive indices for layer {}/{} (sampled at {} points)",
263 layer, solver->lcount, refine * nN);
264 #endif
265
266 if (isnan(lam))
267 throw BadInput(SOLVER->getId(), "no wavelength given: specify 'lam' or 'lam0'");
268
269 if (epsilon_connected && solver->lcomputed[layer]) {
270 SOLVER->writelog(LOG_DEBUG, "Layer {:d} takes some materials parameters from inEpsilon", layer);
271 if (isnan(glam)) glam = lam;
272 }
273 if (gain_connected && solver->lgained[layer]) {
274 SOLVER->writelog(LOG_DEBUG, "Layer {:d} has gain", layer);
275 if (isnan(glam)) glam = lam;
276 }
277
278 double factor = 1. / double(refine);
279
280 double pl = left + SOLVER->pml.size, pr = right - SOLVER->pml.size;
281 Tensor3<dcomplex> epsl, epsr;
282 if (!periodic) {
283 double Tl = 0., Tr = 0., totalw = 0.;
284 for (size_t i = 0, vl = pil * solver->verts->size(), vr = pir * solver->verts->size(); i != mesh->vert()->size(); ++vl, ++vr, ++i) {
285 if (solver->stack[i] == layer) {
286 double w = (i == 0 || i == mesh->vert()->size()-1)? 1e-6 : solver->vbounds->at(i) - solver->vbounds->at(i-1);
287 Tl += w * temperature[vl]; Tr += w * temperature[vr]; totalw += w;
288 }
289 }
290 Tl /= totalw; Tr /= totalw;
291 {
292 OmpLockGuard lock; // this must be declared before `material` to guard its destruction
293 auto material = geometry->getMaterial(vec(pl,maty));
294 lock = material->lock();
295 epsl = geometry->getMaterial(vec(pl,maty))->Eps(lam, Tl);
296 if (isnan(epsl))
297 throw BadInput(solver->getId(), "complex permittivity tensor (Eps) for {} is NaN at lam={}nm and T={}K",
298 material->name(), lam, Tl);
299 }{
300 OmpLockGuard lock; // this must be declared before `material` to guard its destruction
301 auto material = geometry->getMaterial(vec(pr,maty));
302 lock = material->lock();
303 epsr = geometry->getMaterial(vec(pr,maty))->Eps(lam, Tr);
304 if (isnan(epsr))
305 throw BadInput(solver->getId(), "complex permittivity tensor (Eps) for {} is NaN at lam={}nm and T={}K",
306 material->name(), lam, Tr);
307 }
308 }
309
310 // Make space for the result
311 if (polarization == E_LONG) {
312 coeffs[layer].zz.reset(nN, 0.);
313 } else {
314 coeffs[layer].rxx.reset(nN, 0.);
315 coeffs[layer].yy.reset(nN, 0.);
316 }
317
318 // Average material parameters
319 for (size_t i = 0; i != nN; ++i) {
320 for (size_t j = refine*i, end = refine*(i+1); j != end; ++j) {
321 Tensor3<dcomplex> eps = getEpsilon(geometry, layer, maty, lam, glam, j);
322
323 // Add PMLs
324 if (!periodic) {
325 if (j < pil) {
326 double h = (pl - mesh->tran()->at(j)) / SOLVER->pml.size;
327 dcomplex sy(1. + (SOLVER->pml.factor-1.)*pow(h, SOLVER->pml.order));
328 eps = Tensor3<dcomplex>(epsl.c00*sy, epsl.c11/sy, epsl.c22*sy);
329 } else if (j > pir) {
330 double h = (mesh->tran()->at(j) - pr) / SOLVER->pml.size;
331 dcomplex sy(1. + (SOLVER->pml.factor-1.)*pow(h, SOLVER->pml.order));
332 eps = Tensor3<dcomplex>(epsr.c00*sy, epsr.c11/sy, epsr.c22*sy);
333 }
334 }
335
336 if (polarization == E_LONG) {
337 if (eps.c01 != 0.)
338 throw BadInput(solver->getId(), "polarization can be specified only for diagonal complex permittivity tensor (Eps)");
339 coeffs[layer].zz[i] += eps.c00;
340 } else if (polarization == E_TRAN) {
341 if (eps.c01 != 0.)
342 throw BadInput(solver->getId(), "polarization can be specified only for diagonal complex permittivity tensor (Eps)");
343 coeffs[layer].rxx[i] += 1./eps.c11;
344 coeffs[layer].yy[i] += eps.c22;
345 } else {
346 dcomplex rm;
347 bool nd = false;
348 if (eps.c01 != 0. || eps.c10 != 0.) {
349 if (!is_zero(eps.c01 - conj(eps.c10)))
350 throw BadInput(solver->getId(), "complex permittivity tensor (Eps) must be Hermitian for this solver");
351 if (epsilon_diagonal) {
352 if (symmetric())
353 throw BadInput(solver->getId(), "symmetry can be specified only for diagonal complex permittivity tensor (Eps)");
354 coeffs[layer].zx.reset(nN, 0.);
355 epsilon_diagonal = false;
356 }
357 nd = true;
358 rm = 1. / (eps.c00*eps.c11 - eps.c01.real()*eps.c01.real() - eps.c01.imag()*eps.c01.imag());
359 coeffs[layer].zx[i] += eps.c01;
360 };
361 if (eps.c00 != eps.c22 || eps.c11 != eps.c22 || !epsilon_isotropic) {
362 if (epsilon_isotropic) {
363 coeffs[layer].zz = coeffs[layer].yy.copy();
364 epsilon_isotropic = false;
365 }
366 coeffs[layer].zz[i] += eps.c00;
367 }
368 coeffs[layer].rxx[i] += nd? rm*eps.c00 : 1./eps.c11;
369
370 coeffs[layer].yy[i] += eps.c22;
371 }
372 }
373
374 if (polarization == E_LONG) {
375 coeffs[layer].zz[i] *= factor;
376 } else {
377 coeffs[layer].rxx[i] *= factor;
378 coeffs[layer].yy[i] *= factor;
379 if (!epsilon_isotropic) {
380 coeffs[layer].zz[i] *= factor;
381 }
382 if (!epsilon_diagonal) {
383 coeffs[layer].zx[i] *= factor;
384 }
385 }
386 }
387
388 // Check if the layer is uniform
389 if (periodic) {
390 diagonals[layer] = true;
391 if (polarization == E_LONG) {
392 for (size_t i = 1; i != nN; ++i)
393 if (!is_zero(coeffs[layer].zz[i] - coeffs[layer].zz[0])) {
394 diagonals[layer] = false; break;
395 }
396 } else if (polarization == E_TRAN) {
397 if (epsilon_isotropic) {
398 for (size_t i = 1; i != nN; ++i)
399 if (!is_zero(coeffs[layer].yy[i] - coeffs[layer].yy[0])) {
400 diagonals[layer] = false; break;
401 }
402 } else {
403 for (size_t i = 1; i != nN; ++i)
404 if (!(is_zero(coeffs[layer].zz[i] - coeffs[layer].zz[0]) && is_zero(coeffs[layer].rxx[i] - coeffs[layer].rxx[0]) && is_zero(coeffs[layer].yy[i] - coeffs[layer].yy[0]))) {
405 diagonals[layer] = false; break;
406 }
407 }
408 if (!epsilon_diagonal) {
409 for (size_t i = 0; i != nN; ++i)
410 if (!is_zero(coeffs[layer].zx[i])) {
411 diagonals[layer] = false; break;
412 }
413 }
414 } else {
415 if (epsilon_isotropic) {
416 for (size_t i = 1; i != nN; ++i)
417 if (!is_zero(coeffs[layer].yy[i] - coeffs[layer].yy[0])) {
418 diagonals[layer] = false; break;
419 }
420 } else
421 diagonals[layer] = false;
422 }
423 } else {
424 diagonals[layer] = false;
425 }
426
427 if (diagonals[layer]) {
428 SOLVER->writelog(LOG_DETAIL, "Layer {0} is uniform", layer);
429 size_t n1 = nN - 1;
430 if (polarization == E_LONG) {
431 std::fill_n(coeffs[layer].zz.data()+1, n1, 0.);
432 } else {
433 std::fill_n(coeffs[layer].rxx.data()+1, n1, 0.);
434 std::fill_n(coeffs[layer].yy.data()+1, n1, 0.);
435 if (epsilon_isotropic) {
436 if (polarization != E_TRAN) {
437 coeffs[layer].zz = coeffs[layer].yy;
438 }
439 } else {
440 std::fill_n(coeffs[layer].zz.data()+1, n1, 0.);
441 }
442 if (!epsilon_diagonal) {
443 std::fill_n(coeffs[layer].zx.data()+1, n1, 0.);
444 }
445 }
446 } else {
447 // Perform FFT
448 if (polarization == E_LONG) {
449 matFFT.execute(coeffs[layer].zz.data());
450 } else {
451 matFFT.execute(coeffs[layer].rxx.data());
452 matFFT.execute(coeffs[layer].yy.data());
453 if (epsilon_isotropic) {
454 if (polarization != E_TRAN) {
455 coeffs[layer].zz = coeffs[layer].yy;
456 }
457 } else {
458 matFFT.execute(coeffs[layer].zz.data());
459 }
460 if (!epsilon_diagonal) {
461 matFFT.execute(coeffs[layer].zx.data());
462 }
463 }
464 }
465
466 } else {
467 #if defined(OPENMP_FOUND) // && !defined(NDEBUG)
468 SOLVER->writelog(LOG_DETAIL, "Getting refractive indices for layer {}/{} in thread {}",
469 layer, solver->lcount, omp_get_thread_num());
470 #else
471 SOLVER->writelog(LOG_DETAIL, "Getting refractive indices for layer {}/{}", layer, solver->lcount);
472 #endif
473
474 if (isnan(lam))
475 throw BadInput(SOLVER->getId(), "no wavelength given: specify 'lam' or 'lam0'");
476
477 if (epsilon_connected && solver->lcomputed[layer]) {
478 SOLVER->writelog(LOG_DEBUG, "Layer {:d} takes some materials parameters from inEpsilon", layer);
479 if (isnan(glam)) glam = lam;
480 }
481 if (gain_connected && solver->lgained[layer]) {
482 SOLVER->writelog(LOG_DEBUG, "Layer {:d} has gain", layer);
483 if (isnan(glam)) glam = lam;
484 }
485
486 size_t mn = mesh->tran()->size();
487 const Tensor3<dcomplex> eps0 = getEpsilon(geometry, layer, maty, lam, glam, mn-1);
488 bool nd = eps0.c01 != 0.;
489 dcomplex rm;
490 if (nd) {
491 if (polarization != E_UNSPECIFIED)
492 throw BadInput(solver->getId(), "polarization can be specified only for diagonal complex permittivity tensor (Eps)");
493 rm = 1. / (eps0.c00*eps0.c11 - eps0.c01.real()*eps0.c01.real() - eps0.c01.imag()*eps0.c01.imag());
494 }
495 const Tensor3<dcomplex> reps0 = nd? Tensor3<dcomplex>(rm*eps0.c11, rm*eps0.c00, 1./eps0.c22, -rm*eps0.c01) :
496 Tensor3<dcomplex>(1./eps0.c00, 1./eps0.c11, 1./eps0.c22);
497
498 if (polarization == E_LONG) {
499 coeffs[layer].zz.reset(nN, 0.); coeffs[layer].zz[0] = eps0.c00;
500 } else {
501 coeffs[layer].rxx.reset(nN, 0.); coeffs[layer].rxx[0] = reps0.c11;
502 coeffs[layer].yy.reset(nN, 0.); coeffs[layer].yy[0] = eps0.c22;
503 if (polarization != E_TRAN) {
504 if (eps0.c00 != eps0.c22 || eps0.c11 != eps0.c22) {
505 coeffs[layer].zz.reset(nN, 0.); coeffs[layer].zz[0] = eps0.c00;
506 } else {
507 coeffs[layer].zz = coeffs[layer].yy;
508 }
509 }
510 }
511
512 diagonals[layer] = true;
513
514 Tensor3<dcomplex> eps = getEpsilon(geometry, layer, maty, lam, glam, 0), reps;
515
516 double l, r = 0.;
517 const ptrdiff_t di = (mesh->tran()->size() == original_mesh->size()+1)? 1 : 0;
518 const int start = symmetric()? 0 : -int(nN)/2, end = symmetric()? nN : int(nN+1)/2;
519 const double b = 2*PI / L;
520 for (size_t i = 1; i < mn; ++i) {
521 Tensor3<dcomplex> eps1 = getEpsilon(geometry, layer, maty, lam, glam, i);
522 if (!eps1.equals(eps)) {
523 nd = eps.c01 != 0.;
524 if (nd) {
525 if (polarization != E_UNSPECIFIED)
526 throw BadInput(solver->getId(), "polarization can be specified only for diagonal complex permittivity tensor (Eps)");
527 rm = 1. / (eps.c00*eps.c11 - eps.c01.real()*eps.c01.real() - eps.c01.imag()*eps.c01.imag());
528 Tensor3<dcomplex>(rm*eps.c11, rm*eps.c00, 1./eps.c22, -rm*eps.c01);
529 } else
530 reps = Tensor3<dcomplex>(1./eps.c00, 1./eps.c11, 1./eps.c22);
531 diagonals[layer] = false;
532 l = r;
533 r = original_mesh->at(i-di) - left;
534 if (polarization == E_LONG) {
535 add_coeffs(start, end, b, l, r, coeffs[layer].zz, eps.c00 - eps0.c00);
536 } else if (polarization == E_TRAN) {
537 add_coeffs(start, end, b, l, r, coeffs[layer].rxx, 1./eps.c11 - reps0.c11);
538 add_coeffs(start, end, b, l, r, coeffs[layer].yy, eps.c22 - eps0.c22);
539 } else {
540 if (eps.c00 != eps.c22 || eps.c11 != eps.c22 || !epsilon_isotropic) {
541 if (epsilon_isotropic) {
542 coeffs[layer].zz = coeffs[layer].yy.copy();
543 epsilon_isotropic = false;
544 }
545 add_coeffs(start, end, b, l, r, coeffs[layer].zz, eps.c00 - eps0.c00);
546 }
547 add_coeffs(start, end, b, l, r, coeffs[layer].rxx, reps.c11 - reps0.c11);
548 add_coeffs(start, end, b, l, r, coeffs[layer].yy, eps.c22 - eps0.c22);
549 if (eps.c01 != 0.) {
550 if (epsilon_diagonal) {
551 if (symmetric())
552 throw BadInput(solver->getId(), "symmetry can be specified only for diagonal complex permittivity tensor (Eps)");
553 coeffs[layer].zx.reset(nN, 0.);
554 epsilon_diagonal = false;
555 }
556 add_coeffs(start, end, b, l, r, coeffs[layer].zx, eps.c01 - eps0.c01);
557 }
558 }
559 eps = eps1;
560 }
561 }
562 //TODO Add PMLs
563 }
564 // Smooth coefficients
565 if (!diagonals[layer] && SOLVER->smooth) {
566 double bb4 = PI / L; bb4 *= bb4; // (2π/L)² / 4
567 for (size_t i = 0; i != nN; ++i) {
568 int k = int(i); if (!symmetric() && k > int(nN/2)) k -= int(nN);
569 double s = exp(-SOLVER->smooth * bb4 * k * k);
570 coeffs[layer].yy[i] *= s;
571 if (polarization == E_LONG) {
572 if (!epsilon_isotropic) coeffs[layer].zz[i] *= s;
573 } else {
574 coeffs[layer].rxx[i] *= s;
575 if (!epsilon_isotropic) {
576 coeffs[layer].zz[i] *= s;
577 }
578 if (!epsilon_diagonal) {
579 coeffs[layer].zx[i] *= s;
580 }
581 }
582 }
583 }
584
585 //TODO there is no need to do all the above if only symmetry changes (just do the code below)
586
587 // Compute necessary inverses of Toeplitz matrices
588 TempMatrix temp = getTempMatrix();
589 cmatrix work(N, N, temp.data());
590
591 const int order = int(SOLVER->getSize());
592
593 if (polarization != E_LONG) {
594 if (symmetric()) {
595 // Full symmetric()
596 const bool sel = symmetry == E_LONG;
597
598 for (int i = 0; i <= order; ++i)
599 work(i,0) = epsyy(layer,i);
600 for (int j = 1; j <= order; ++j)
601 for (int i = 0; i <= order; ++i)
602 work(i,j) = sel? epsyy(layer,abs(i-j)) + epsyy(layer,i+j) : epsyy(layer,abs(i-j)) - epsyy(layer,i+j);
603 coeff_matrices[layer].reyy.reset(N, N);
604 make_unit_matrix(coeff_matrices[layer].reyy);
605 invmult(work, coeff_matrices[layer].reyy);
606
607 for (int i = 0; i <= order; ++i)
608 work(i,0) = repsxx(layer,i);
609 for (int j = 1; j <= order; ++j)
610 for (int i = 0; i <= order; ++i)
611 work(i,j) = sel? repsxx(layer,abs(i-j)) - repsxx(layer,i+j) : repsxx(layer,abs(i-j)) + repsxx(layer,i+j);
612 coeff_matrices[layer].exx.reset(N, N);
613 make_unit_matrix(coeff_matrices[layer].exx);
614 invmult(work, coeff_matrices[layer].exx);
615 } else {
616 for (int j = -order; j <= order; ++j) {
617 const size_t jt = iEH(j);
618 for (int i = -order; i <= order; ++i) {
619 const size_t it = iEH(i);
620 work(it,jt) = epsyy(layer,i-j);
621 }
622 }
623 coeff_matrices[layer].reyy.reset(N, N);
624 make_unit_matrix(coeff_matrices[layer].reyy);
625 invmult(work, coeff_matrices[layer].reyy);
626
627 for (int j = -order; j <= order; ++j) {
628 const size_t jt = iEH(j);
629 for (int i = -order; i <= order; ++i) {
630 const size_t it = iEH(i);
631 work(it,jt) = repsxx(layer,i-j);
632 }
633 }
634 coeff_matrices[layer].exx.reset(N, N);
635 make_unit_matrix(coeff_matrices[layer].exx);
636 invmult(work, coeff_matrices[layer].exx);
637 }
638 }
639}
640
641
642LazyData<Tensor3<dcomplex>> ExpansionPW2D::getMaterialEps(size_t l, const shared_ptr<const LevelsAdapter::Level> &level, InterpolationMethod interp)
643{
644 assert(dynamic_pointer_cast<const MeshD<2>>(level->mesh()));
645 auto dest_mesh = static_pointer_cast<const MeshD<2>>(level->mesh());
646 if (interp == INTERPOLATION_DEFAULT || interp == INTERPOLATION_FOURIER) {
647 if (!symmetric()) {
648 return LazyData<Tensor3<dcomplex>>(dest_mesh->size(), [this,l,dest_mesh](size_t i)->Tensor3<dcomplex>{
649 Tensor3<dcomplex> eps(0.);
650 for (int k = -int(nN)/2, end = int(nN+1)/2; k != end; ++k) {
651 size_t j = (k>=0)? k : k + nN;
652 dcomplex ff = exp(2*PI * k * I * (dest_mesh->at(i).c0-left) / (right-left));
653 switch (polarization) {
654 case E_LONG:
655 eps.c00 += coeffs[l].zz[j] * ff;
656 break;
657 case E_TRAN:
658 eps.c11 += coeffs[l].rxx[j] * ff;
659 eps.c22 += coeffs[l].yy[j] * ff;
660 break;
661 case E_UNSPECIFIED:
662 eps.c00 += coeffs[l].zz[j] * ff;
663 eps.c11 += coeffs[l].rxx[j] * ff;
664 eps.c22 += coeffs[l].yy[j] * ff;
665 if (coeffs[l].zx) eps.c01 += coeffs[l].zx[k] * ff;
666 break;
667 }
668 }
669 switch (polarization) {
670 case E_LONG:
671 eps.c22 = eps.c11 = eps.c00;
672 break;
673 case E_TRAN:
674 eps.c00 = eps.c22;
675 case E_UNSPECIFIED:
676 eps.c11 = 1. / eps.c11;
677 }
678 return eps;
679 });
680 } else {
681 return LazyData<Tensor3<dcomplex>>(dest_mesh->size(), [this,l,dest_mesh](size_t i)->Tensor3<dcomplex>{
682 Tensor3<dcomplex> eps(0.);
683 for (std::size_t k = 0; k != nN; ++k) {
684 dcomplex ff = (k? 2. : 1.) * cos(PI * double(k) * dest_mesh->at(i).c0 / (right-left));
685 switch (polarization) {
686 case E_LONG:
687 eps.c00 += coeffs[l].zz[k] * ff;
688 break;
689 case E_TRAN:
690 eps.c11 += coeffs[l].rxx[k] * ff;
691 eps.c22 += coeffs[l].yy[k] * ff;
692 break;
693 case E_UNSPECIFIED:
694 eps.c00 += coeffs[l].zz[k] * ff;
695 eps.c11 += coeffs[l].rxx[k] * ff;
696 eps.c22 += coeffs[l].yy[k] * ff;
697 break;
698 }
699 }
700 switch (polarization) {
701 case E_LONG:
702 eps.c22 = eps.c11 = eps.c00;
703 break;
704 case E_TRAN:
705 eps.c00 = eps.c22;
706 case E_UNSPECIFIED:
707 eps.c11 = 1. / eps.c11;
708 }
709 return eps;
710 });
711 }
712 } else {
713 DataVector<Tensor3<dcomplex>> params(symmetric()? nN : nN+1, Tensor3<dcomplex>(0.));
714 FFT::Backward1D fft(9, int(nN), symmetric()? SOLVER->dct2()? FFT::SYMMETRY_EVEN_2 : FFT::SYMMETRY_EVEN_1 : FFT::SYMMETRY_NONE);
715 if (symmetry == E_LONG) {
716 for (size_t i = 0; i != nN; ++i) params[i].c00 = coeffs[l].zz[i];
717 fft.execute(reinterpret_cast<dcomplex*>(params.data(), 1));
718 for (Tensor3<dcomplex>& eps: params) {
719 eps.c22 = eps.c11 = eps.c00;
720 }
721 } else {
722 for (size_t i = 0; i != nN; ++i) {
723 params[i].c11 = coeffs[l].rxx[i];
724 params[i].c22 = coeffs[l].yy[i];
725 }
726 fft.execute(reinterpret_cast<dcomplex*>(params.data())+1, 1);
727 fft.execute(reinterpret_cast<dcomplex*>(params.data())+2, 1);
728 if (coeffs[l].zx) {
729 for (size_t i = 0; i != nN; ++i) params[i].c01 = coeffs[l].zx[i];
730 fft.execute(reinterpret_cast<dcomplex*>(params.data())+3, 1);
731 for (size_t i = 0; i != nN; ++i) params[i].c10 = conj(params[i].c01);
732 }
733 if (symmetry == E_TRAN || coeffs[l].zz.data() == coeffs[l].yy.data()) {
734 for (Tensor3<dcomplex>& eps: params) {
735 eps.c00 = eps.c22;
736 eps.c11 = 1. / eps.c11;
737 }
738 } else {
739 for (size_t i = 0; i != nN; ++i) params[i].c00 = coeffs[l].zz[i];
740 fft.execute(reinterpret_cast<dcomplex*>(params.data()), 1);
741 for (Tensor3<dcomplex>& eps: params) {
742 eps.c11 = 1. / eps.c11;
743 }
744 }
745 }
746 shared_ptr<RegularAxis> cmesh = plask::make_shared<RegularAxis>();
747 if (symmetric()) {
748 if (SOLVER->dct2()) {
749 double dx = 0.5 * right / double(nN);
750 cmesh->reset(dx, right-dx, nN);
751 } else {
752 cmesh->reset(0., right, nN);
753 }
754 } else {
755 cmesh->reset(left, right, nN+1);
756 params[nN] = params[0];
757 }
758 auto src_mesh = plask::make_shared<RectangularMesh<2>>(cmesh, plask::make_shared<RegularAxis>(level->vpos(), level->vpos(), 1));
759 return interpolate(src_mesh, params, dest_mesh, interp,
760 InterpolationFlags(SOLVER->getGeometry(),
761 symmetric()? InterpolationFlags::Symmetry::POSITIVE : InterpolationFlags::Symmetry::NO,
762 InterpolationFlags::Symmetry::NO)
763 );
764 }
765}
766
767void ExpansionPW2D::make_permeability_matrices(cmatrix& work)
768{
769 coeff_matrix_rmyy.reset(N, N);
770
771 int order = int(SOLVER->getSize());
772
773 if (symmetric()) {
774 const bool sym = symmetry == E_LONG;
775 for (int i = 0; i <= order; ++i)
776 work(i,0) = muyy(i);
777 for (int j = 1; j <= order; ++j)
778 for (int i = 0; i <= order; ++i)
779 work(i,j) = sym? muyy(abs(i-j)) - muyy(i+j) : muyy(abs(i-j)) + muyy(i+j);
780 make_unit_matrix(coeff_matrix_rmyy);
781 invmult(work, coeff_matrix_rmyy);
782
783 if (polarization != E_TRAN) {
784 coeff_matrix_mxx.reset(N, N);
785 for (int i = 0; i <= order; ++i)
786 work(i,0) = rmuxx(i);
787 for (int j = 1; j <= order; ++j)
788 for (int i = 0; i <= order; ++i)
789 work(i,j) = sym? rmuxx(abs(i-j)) + rmuxx(i+j) : rmuxx(abs(i-j)) - rmuxx(i+j);
790 make_unit_matrix(coeff_matrix_mxx);
791 invmult(work, coeff_matrix_mxx);
792 }
793 } else {
794 for (int j = -order; j <= order; ++j) {
795 const size_t jt = iEH(j);
796 for (int i = -order; i <= order; ++i) {
797 const size_t it = iEH(i);
798 work(it,jt) = muyy(i-j);
799 }
800 }
801 make_unit_matrix(coeff_matrix_rmyy);
802 invmult(work, coeff_matrix_rmyy);
803
804 if (polarization != E_TRAN) {
805 coeff_matrix_mxx.reset(N, N);
806 for (int j = -order; j <= order; ++j) {
807 const size_t jt = iEH(j);
808 for (int i = -order; i <= order; ++i) {
809 const size_t it = iEH(i);
810 work(it,jt) = rmuxx(i-j);
811 }
812 }
813 make_unit_matrix(coeff_matrix_mxx);
814 invmult(work, coeff_matrix_mxx);
815 }
816 }
817}
818
819
820void ExpansionPW2D::getMatrices(size_t l, cmatrix& RE, cmatrix& RH)
821{
822 assert(initialized);
823 if (isnan(k0)) throw BadInput(SOLVER->getId(), "wavelength or k0 not set");
824 if (isinf(k0.real())) throw BadInput(SOLVER->getId(), "wavelength must not be 0");
825
826 if (coeff_matrices.empty()) coeff_matrices.resize(solver->lcount);
827
828 dcomplex beta{ this->beta.real(), this->beta.imag() - SOLVER->getMirrorLosses(this->beta.real()/k0.real()) };
829
830 const int order = int(SOLVER->getSize());
831 dcomplex rk0 = 1. / k0, k02 = k0*k0;
832 double b = 2.*PI / (right-left) * (symmetric()? 0.5 : 1.0);
833 double bb = b * b;
834
835 // Ez represents -Ez
836
837 if (separated()) {
838 if (symmetric()) {
839 // Separated symmetric()
840 if (polarization == E_LONG) { // Ez & Hx
841 const bool sym = symmetry == E_LONG;
842 if (!periodic) {
843 coeff_matrix_rmyy.copyto(RE);
844 coeff_matrix_mxx.copyto(RH);
845 } else {
846 make_unit_matrix(RE);
847 make_unit_matrix(RH);
848 }
849 for (int i = 0; i <= order; ++i) {
850 dcomplex gi = - rk0 * bb * double(i);
851 RE(i,0) = k0 * epszz(l,i);
852 RH(i,0) *= k0;
853 for (int j = 1; j <= order; ++j) {
854 RE(i,j) = gi * double(j) * RE(i,j) + k0 * (sym? epszz(l,abs(i-j)) + epszz(l,i+j) : epszz(l,abs(i-j)) - epszz(l,i+j));
855 RH(i,j) *= k0;
856 }
857 // Ugly hack to avoid singularity
858 if (RE(i,i) == 0.) RE(i,i) = 1e-32;
859 if (RH(i,i) == 0.) RH(i,i) = 1e-32;
860 }
861 } else { // Ex & Hz
862 const bool sym = symmetry == E_TRAN;
863 coeff_matrices[l].exx.copyto(RE);
864 coeff_matrices[l].reyy.copyto(RH);
865 for (int i = 0; i <= order; ++i) {
866 const dcomplex gi = - rk0 * bb * double(i);
867 RE(i,0) *= k0;
868 RH(i,0) = k0 * muzz(i);
869 for (int j = 1; j <= order; ++j) {
870 RE(i,j) *= k0;
871 RH(i,j) = gi * double(j) * RH(i,j) + k0 * (sym? muzz(abs(i-j)) + muzz(i+j) : muzz(abs(i-j)) - muzz(i+j));
872 }
873 // Ugly hack to avoid singularity
874 if (RE(i,i) == 0.) RE(i,i) = 1e-32;
875 if (RH(i,i) == 0.) RH(i,i) = 1e-32;
876 }
877 }
878 } else {
879 // Separated asymmetric()
880 if (polarization == E_LONG) { // Ez & Hx
881 if (!periodic) {
882 coeff_matrix_rmyy.copyto(RE);
883 coeff_matrix_mxx.copyto(RH);
884 } else {
885 make_unit_matrix(RE);
886 make_unit_matrix(RH);
887 }
888 for (int i = -order; i <= order; ++i) {
889 const dcomplex gi = - rk0 * (b * double(i) - ktran);
890 const size_t it = iEH(i);
891 for (int j = -order; j <= order; ++j) {
892 const size_t jt = iEH(j);
893 RE(it,jt) = gi * (b * double(j) - ktran) * RE(it,jt) + k0 * epszz(l,i-j);
894 RH(it,jt) *= k0;
895 }
896 // Ugly hack to avoid singularity
897 if (RE(it,it) == 0.) RE(it,it) = 1e-32;
898 if (RH(it,it) == 0.) RH(it,it) = 1e-32;
899 }
900 } else { // Ex & Hz
901 coeff_matrices[l].exx.copyto(RE);
902 coeff_matrices[l].reyy.copyto(RH);
903 for (int i = -order; i <= order; ++i) {
904 const dcomplex gi = - rk0 * (b * double(i) - ktran);
905 const size_t it = iEH(i);
906 for (int j = -order; j <= order; ++j) {
907 const size_t jt = iEH(j);
908 RE(it,jt) *= k0;
909 RH(it,jt) = gi * (b * double(j) - ktran) * RH(it,jt) + k0 * muzz(i-j);
910 }
911 // Ugly hack to avoid singularity
912 if (RE(it,it) == 0.) RE(it,it) = 1e-32;
913 if (RH(it,it) == 0.) RH(it,it) = 1e-32;
914 }
915 }
916 }
917 } else {
918 // work matrix is 2N×2N, so we can use its space for four N×N matrices
919 const size_t NN = N*N;
920 TempMatrix temp = getTempMatrix();
921 cmatrix work(temp);
922 cmatrix workxx, workyy;
923 if (symmetric()) {
924 // Full symmetric()
925 const bool sel = symmetry == E_LONG;
926 workyy = coeff_matrices[l].reyy;
927 if (!periodic) {
928 workxx = coeff_matrix_mxx;
929 } else {
930 workxx = cmatrix(N, N, work.data()+NN);
931 make_unit_matrix(workxx);
932 }
933 for (int i = 0; i <= order; ++i) {
934 const dcomplex gi = b * double(i) - ktran;
935 const size_t iex = iEx(i), iez = iEz(i);
936 for (int j = 0; j <= order; ++j) {
937 int ijp = abs(i-j), ijn = i+j;
938 dcomplex gj = b * double(j) - ktran;
939 const size_t jhx = iHx(j), jhz = iHz(j);
940 RH(iex,jhz) = - rk0 * gi * gj * workyy(i,j) +
941 k0 * (j == 0? muzz(i) : sel? muzz(ijp) - muzz(ijn) : muzz(ijp) + muzz(ijn));
942 RH(iex,jhx) = - rk0 * beta* gi * workyy(i,j);
943 RH(iez,jhz) = - rk0 * beta* gj * workyy(i,j);
944 RH(iez,jhx) = - rk0 * beta*beta * workyy(i,j) + k0 * workxx(i,j);
945 }
946 // Ugly hack to avoid singularity
947 if (RH(iex,iex) == 0.) RH(iex,iex) = 1e-32;
948 if (RH(iez,iez) == 0.) RH(iez,iez) = 1e-32;
949 }
950 workxx = coeff_matrices[l].exx;
951 if (!periodic) {
952 workyy = coeff_matrix_rmyy;
953 } else {
954 workyy = cmatrix(N, N, work.data()+2*NN);
955 make_unit_matrix(workyy);
956 }
957 for (int i = 0; i <= order; ++i) {
958 const dcomplex gi = b * double(i) - ktran;
959 const size_t ihx = iHx(i), ihz = iHz(i);
960 for (int j = 0; j <= order; ++j) {
961 int ijp = abs(i-j), ijn = i+j;
962 dcomplex gj = b * double(j) - ktran;
963 const size_t jex = iEx(j), jez = iEz(j);
964 RE(ihz,jex) = - rk0 * beta*beta * workyy(i,j) + k0 * workxx(i,j);
965 RE(ihz,jez) = rk0 * beta* gj * workyy(i,j);
966 RE(ihx,jex) = rk0 * beta* gi * workyy(i,j);
967 RE(ihx,jez) = - rk0 * gi * gj * workyy(i,j) +
968 k0 * (j == 0? epszz(l,i) : sel? epszz(l,ijp) + epszz(l,ijn) : epszz(l,ijp) - epszz(l,ijn));
969 }
970 // Ugly hack to avoid singularity
971 if (RE(ihx,ihx) == 0.) RE(ihx,ihx) = 1e-32;
972 if (RE(ihz,ihz) == 0.) RE(ihz,ihz) = 1e-32;
973 }
974 } else {
975 // Full asymmetric()
976 workyy = coeff_matrices[l].reyy;
977 if (!periodic) {
978 workxx = coeff_matrix_mxx;
979 } else {
980 workxx = cmatrix(N, N, work.data()+NN);
981 make_unit_matrix(workxx);
982 }
983 for (int i = -order; i <= order; ++i) {
984 const dcomplex gi = b * double(i) - ktran;
985 const size_t iex = iEx(i), iez = iEz(i), it = iEH(i);
986 for (int j = -order; j <= order; ++j) {
987 int ij = i-j; dcomplex gj = b * double(j) - ktran;
988 const size_t jhx = iHx(j), jhz = iHz(j), jt = iEH(j);
989 RH(iex,jhz) = - rk0 * gi * gj * workyy(it,jt) + k0 * muzz(ij);
990 RH(iex,jhx) = - rk0 * beta* gi * workyy(it,jt);
991 RH(iez,jhz) = - rk0 * beta* gj * workyy(it,jt);
992 RH(iez,jhx) = - rk0 * beta*beta * workyy(it,jt) + k0 * workxx(it,jt);
993 }
994 // Ugly hack to avoid singularity
995 if (RH(iex,iex) == 0.) RH(iex,iex) = 1e-32;
996 if (RH(iez,iez) == 0.) RH(iez,iez) = 1e-32;
997 }
998 workxx = coeff_matrices[l].exx;
999 if (!periodic) {
1000 workyy = coeff_matrix_rmyy;
1001 } else {
1002 workyy = cmatrix(N, N, work.data()+2*NN);
1003 make_unit_matrix(workyy);
1004 }
1005 for (int i = -order; i <= order; ++i) {
1006 const dcomplex gi = b * double(i) - ktran;
1007 const size_t ihx = iHx(i), ihz = iHz(i), it = iEH(i);
1008 for (int j = -order; j <= order; ++j) {
1009 int ij = i-j; dcomplex gj = b * double(j) - ktran;
1010 const size_t jex = iEx(j), jez = iEz(j), jt = iEH(j);
1011 RE(ihz,jex) = - rk0 * beta*beta * workyy(it,jt) + k0 * workxx(it,jt);
1012 RE(ihz,jez) = rk0 * beta* gj * workyy(it,jt);
1013 RE(ihx,jex) = rk0 * beta* gi * workyy(it,jt);
1014 RE(ihx,jez) = - rk0 * gi * gj * workyy(it,jt) + k0 * epszz(l,ij);
1015 if (epszx(l)) {
1016 RE(ihx,jex) -= k0 * epszx(l,ij);
1017 RE(ihz,jez) -= k0 * epsxz(l,ij);
1018 }
1019 }
1020 // Ugly hack to avoid singularity
1021 if (RE(ihx,ihx) == 0.) RE(ihx,ihx) = 1e-32;
1022 if (RE(ihz,ihz) == 0.) RE(ihz,ihz) = 1e-32;
1023 }
1024 }
1025 }
1026}
1027
1028
1029void ExpansionPW2D::prepareField()
1030{
1031 if (field_interpolation == INTERPOLATION_DEFAULT) field_interpolation = INTERPOLATION_FOURIER;
1032 if (symmetric()) {
1033 field.reset(N);
1034 if (field_interpolation != INTERPOLATION_FOURIER) {
1035 Component sym = (which_field == FIELD_E || !symmetry)? symmetry : Component(3-symmetry);
1036 int df = SOLVER->dct2()? 0 : 4;
1037 fft_x = FFT::Backward1D(3, N, FFT::Symmetry(sym+df)); // tran
1038 fft_yz = FFT::Backward1D(3, N, FFT::Symmetry(3-sym+df)); // long
1039 }
1040 } else {
1041 field.reset(N + 1);
1042 if (field_interpolation != INTERPOLATION_FOURIER)
1043 fft_x = FFT::Backward1D(3, N, FFT::SYMMETRY_NONE);
1044 }
1045}
1046
1047void ExpansionPW2D::cleanupField()
1048{
1049 field.reset();
1050 fft_x = FFT::Backward1D();
1051 fft_yz = FFT::Backward1D();
1052}
1053
1054// TODO fields must be carefully verified
1055
1056LazyData<Vec<3,dcomplex>> ExpansionPW2D::getField(size_t l, const shared_ptr<const typename LevelsAdapter::Level> &level, const cvector& E, const cvector& H)
1057{
1058 Component sym = (which_field == FIELD_E)? symmetry : Component((3-symmetry) % 3);
1059
1060 dcomplex beta{ this->beta.real(), this->beta.imag() - SOLVER->getMirrorLosses(this->beta.real()/k0.real()) };
1061
1062 const int order = int(SOLVER->getSize());
1063 double b = 2.*PI / (right-left) * (symmetric()? 0.5 : 1.0);
1064 assert(dynamic_pointer_cast<const MeshD<2>>(level->mesh()));
1065 auto dest_mesh = static_pointer_cast<const MeshD<2>>(level->mesh());
1066 double vpos = level->vpos();
1067
1068 int dx = (symmetric() && field_interpolation != INTERPOLATION_FOURIER && sym != E_TRAN)? 1 : 0; // 1 for sin expansion of tran component
1069 int dz = (symmetric() && field_interpolation != INTERPOLATION_FOURIER && sym != E_LONG)? 1 : 0; // 1 for sin expansion of long component
1070
1071 TempMatrix temp = getTempMatrix();
1072 cvector work(temp.data(), N);
1073
1074 if (which_field == FIELD_E) {
1075 if (polarization == E_LONG) {
1076 for (int i = symmetric()? 0 : -order; i <= order; ++i) {
1077 size_t ieh = iEH(i);
1078 field[ieh].tran() = field[ieh].vert() = 0.;
1079 if (ieh != 0 || !dz) field[ieh-dz].lon() = - E[ieh];
1080 }
1081 } else if (polarization == E_TRAN) {
1082 for (int i = symmetric()? 0 : -order; i <= order; ++i) {
1083 size_t ieh = iEH(i);
1084 field[ieh].lon() = 0.;
1085 if (ieh != 0 || !dx)
1086 field[ieh-dx].tran() = E[ieh];
1087 if (ieh != 0 || !dz) {
1088 field[ieh-dz].vert() = 0.;
1089 for (int j = symmetric()? 0 : -order; j <= order; ++j) {
1090 size_t jeh = iEH(j);
1091 field[ieh-dz].vert() += coeff_matrices[l].reyy(ieh,jeh) * (b*double(j)-ktran) * H[jeh];
1092 }
1093 field[ieh-dz].vert() /= k0;
1094 }
1095 }
1096 } else {
1097 for (int i = symmetric()? 0 : -order; i <= order; ++i) {
1098 size_t ieh = iEH(i);
1099 if (ieh != 0 || !dx)
1100 field[ieh-dx].tran() = E[iEx(i)];
1101 if (ieh != 0 || !dz) {
1102 field[ieh-dz].lon() = - E[iEz(i)];
1103 field[ieh-dz].vert() = 0.;
1104 for (int j = symmetric()? 0 : -order; j <= order; ++j) {
1105 field[ieh-dz].vert() += coeff_matrices[l].reyy(ieh,iEH(j))
1106 * ((b*double(j)-ktran) * H[iHz(j)] - beta * H[iHx(j)]);
1107 }
1108 field[ieh-dz].vert() /= k0;
1109 }
1110 }
1111 }
1112 } else { // which_field == FIELD_H
1113 if (polarization == E_TRAN) { // polarization == H_LONG
1114 for (int i = symmetric()? 0 : -order; i <= order; ++i) {
1115 size_t ieh = iEH(i);
1116 field[ieh].tran() = field[ieh].vert() = 0.;
1117 if (ieh != 0 || !dz) field[ieh-dz].lon() = H[ieh];
1118 }
1119 } else if (polarization == E_LONG) { // polarization == H_TRAN
1120 for (int i = symmetric()? 0 : -order; i <= order; ++i) {
1121 size_t ieh = iEH(i);
1122 field[ieh].lon() = 0.;
1123 if (ieh != 0 || !dx)
1124 field[ieh-dx].tran() = H[ieh];
1125 if (ieh != 0 || !dz) {
1126 if (periodic)
1127 field[ieh-dz].vert() = - (b*double(i)-ktran) * E[ieh];
1128 else {
1129 field[ieh-dz].vert() = 0.;
1130 for (int j = symmetric()? 0 : -order; j <= order; ++j) {
1131 size_t jeh = iEH(j);
1132 field[ieh-dz].vert() -= coeff_matrix_rmyy(ieh,jeh) * (b*double(j)-ktran) * E[jeh];
1133 }
1134 }
1135 field[ieh-dz].vert() /= k0;
1136 }
1137 }
1138 } else {
1139 for (int i = symmetric()? 0 : -order; i <= order; ++i) {
1140 size_t ieh = iEH(i);
1141 if (ieh != 0 || !dx)
1142 field[ieh-dx].tran() = H[iHx(i)];
1143 if (ieh != 0 || !dz) {
1144 field[ieh-dz].lon() = H[iHz(i)];
1145 if (periodic)
1146 field[ieh-dz].vert() = (beta * E[iEx(i)] - (b*double(i)-ktran) * E[iEz(i)]);
1147 else {
1148 field[ieh-dz].vert() = 0.;
1149 for (int j = symmetric()? 0 : -order; j <= order; ++j)
1150 field[ieh-dz].vert() += coeff_matrix_rmyy(ieh,iEH(j))
1151 * (beta * E[iEx(j)] - (b*double(j)-ktran) * E[iEz(j)]);
1152 }
1153 field[ieh-dz].vert() /= k0;
1154 }
1155 }
1156 }
1157 }
1158
1159 if (dx) { field[field.size()-1].tran() = 0.; }
1160 if (dz) { field[field.size()-1].lon() = 0.; field[field.size()-1].vert() = 0.; }
1161
1162 if (field_interpolation == INTERPOLATION_FOURIER) {
1163 DataVector<Vec<3,dcomplex>> result(dest_mesh->size());
1164 double L = right - left;
1165 if (!symmetric()) {
1166 dcomplex B = 2*PI * I / L;
1167 dcomplex ikx = I * ktran;
1168 result.reset(dest_mesh->size(), Vec<3,dcomplex>(0.,0.,0.));
1169 for (int k = -order; k <= order; ++k) {
1170 size_t j = (k>=0)? k : k + N;
1171 dcomplex G = B * double(k) - ikx;
1172 for (size_t i = 0; i != dest_mesh->size(); ++i) {
1173 double x = dest_mesh->at(i)[0];
1174 if (!periodic) x = clamp(x, left, right);
1175 result[i] += field[j] * exp(G * (x-left));
1176 }
1177 }
1178 } else {
1179 double B = PI / L;
1180 result.reset(dest_mesh->size());
1181 for (size_t i = 0; i != dest_mesh->size(); ++i) {
1182 result[i] = field[0];
1183 for (int k = 1; k <= order; ++k) {
1184 double x = dest_mesh->at(i)[0];
1185 if (!periodic) x = clamp(x, -right, right);
1186 double cs = 2. * cos(B * k * x);
1187 dcomplex sn = 2. * I * sin(B * k * x);
1188 if (sym == E_TRAN) {
1189 result[i].lon() += field[k].lon() * sn;
1190 result[i].tran() += field[k].tran() * cs;
1191 result[i].vert() += field[k].vert() * sn;
1192 } else {
1193 result[i].lon() += field[k].lon() * cs;
1194 result[i].tran() += field[k].tran() * sn;
1195 result[i].vert() += field[k].vert() * cs;
1196 }
1197 }
1198 }
1199 }
1200 return result;
1201 } else {
1202 if (symmetric()) {
1203 fft_x.execute(&(field.data()->tran()), 1);
1204 fft_yz.execute(&(field.data()->lon()), 1);
1205 fft_yz.execute(&(field.data()->vert()), 1);
1206 if (sym == E_TRAN) {
1207 for (Vec<3,dcomplex>& f: field) {
1208 f.c0 = dcomplex(-f.c0.imag(), f.c0.real());
1209 f.c2 = dcomplex(-f.c2.imag(), f.c2.real());
1210
1211 }
1212 } else {
1213 for (Vec<3,dcomplex>& f: field) {
1214 f.c1 = dcomplex(-f.c1.imag(), f.c1.real());
1215
1216 }
1217 }
1218 double dx = 0.5 * (right-left) / double(N);
1219 auto src_mesh = plask::make_shared<RectangularMesh<2>>(plask::make_shared<RegularAxis>(left+dx, right-dx, field.size()), plask::make_shared<RegularAxis>(vpos, vpos, 1));
1220 return interpolate(src_mesh, field, dest_mesh, field_interpolation,
1221 InterpolationFlags(SOLVER->getGeometry(),
1222 (sym == E_TRAN)? InterpolationFlags::Symmetry::NPN : InterpolationFlags::Symmetry::PNP,
1223 InterpolationFlags::Symmetry::NO),
1224 false);
1225 } else {
1226 fft_x.execute(reinterpret_cast<dcomplex*>(field.data()));
1227 field[N] = field[0];
1228 auto src_mesh = plask::make_shared<RectangularMesh<2>>(plask::make_shared<RegularAxis>(left, right, field.size()), plask::make_shared<RegularAxis>(vpos, vpos, 1));
1229 LazyData<Vec<3,dcomplex>> interpolated =
1230 interpolate(src_mesh, field, dest_mesh, field_interpolation,
1231 InterpolationFlags(SOLVER->getGeometry(),
1232 InterpolationFlags::Symmetry::NO, InterpolationFlags::Symmetry::NO),
1233 false);
1234 dcomplex ikx = I * ktran;
1235 return LazyData<Vec<3,dcomplex>>(interpolated.size(), [interpolated, dest_mesh, ikx] (size_t i) {
1236 return interpolated[i] * exp(-ikx * dest_mesh->at(i).c0);
1237 });
1238 }
1239 }
1240}
1241
1242
1243double ExpansionPW2D::integratePoyntingVert(const cvector& E, const cvector& H)
1244{
1245 double P = 0.;
1246
1247 const int ord = int(SOLVER->getSize());
1248
1249 if (separated()) {
1250 if (symmetric()) {
1251 for (int i = 0; i <= ord; ++i) {
1252 P += real(E[iEH(i)] * conj(H[iEH(i)]));
1253 }
1254 P = 2. * P - real(E[iEH(0)] * conj(H[iEH(0)]));
1255 } else {
1256 for (int i = -ord; i <= ord; ++i) {
1257 P += real(E[iEH(i)] * conj(H[iEH(i)]));
1258 }
1259 }
1260 } else {
1261 if (symmetric()) {
1262 for (int i = 0; i <= ord; ++i) {
1263 P -= real(E[iEz(i)] * conj(H[iHx(i)]) + E[iEx(i)] * conj(H[iHz(i)]));
1264 }
1265 P = 2. * P + real(E[iEz(0)] * conj(H[iHx(0)]) + E[iEx(0)] * conj(H[iHz(0)]));
1266 } else {
1267 for (int i = -ord; i <= ord; ++i) {
1268 P -= real(E[iEz(i)] * conj(H[iHx(i)]) + E[iEx(i)] * conj(H[iHz(i)]));
1269 }
1270 }
1271 }
1272
1273 double L = SOLVER->geometry->getExtrusion()->getLength();
1274 if (!isinf(L)) P *= L * 1e-6; // µm -> m
1275
1276 return P * (symmetric()? 2 * right : right - left) * 1e-6; // µm -> m
1277}
1278
1279
1280void ExpansionPW2D::getDiagonalEigenvectors(cmatrix& Te, cmatrix& Te1, const cmatrix& RE, const cdiagonal& gamma)
1281{
1282 size_t nr = Te.rows(), nc = Te.cols();
1283 std::fill_n(Te.data(), nr*nc, 0.);
1284 std::fill_n(Te1.data(), nr*nc, 0.);
1285
1286 if (separated()) {
1287 for (std::size_t i = 0; i < nc; i++) {
1288 Te(i,i) = Te1(i,i) = 1.;
1289 }
1290 } else {
1291 // Ensure that for the same gamma E*H [2x2] is diagonal
1292 assert(nc % 2 == 0);
1293 size_t n = nc / 2;
1294 for (std::size_t i = 0; i < n; i++) {
1295 // Compute Te1 = sqrt(RE)
1296 // https://en.wikipedia.org/wiki/Square_root_of_a_2_by_2_matrix
1297 // but after this normalize columns to 1
1298 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);
1299 dcomplex s = sqrt(a*d - b*c);
1300 a += s; d += s;
1301 // Normalize
1302 s = 1. / sqrt(a*a + b*b); a *= s; b *= s;
1303 s = 1. / sqrt(c*c + d*d); c *= s; d *= s;
1304 Te1(2*i, 2*i) = a; Te1(2*i, 2*i+1) = b;
1305 Te1(2*i+1, 2*i) = c; Te1(2*i+1, 2*i+1) = d;
1306 // Invert Te1
1307 s = 1. / (a*d - b*c);
1308 Te(2*i, 2*i) = s * d; Te(2*i, 2*i+1) = - s * b;
1309 Te(2*i+1, 2*i) = - s * c; Te(2*i+1, 2*i+1) = s * a;
1310 }
1311 }
1312}
1313
1314double ExpansionPW2D::integrateField(WhichField field, size_t l, const cmatrix& TE, const cmatrix& TH,
1315 const std::function<std::pair<dcomplex,dcomplex>(size_t, size_t)>& vertical)
1316{
1317 const int order = int(SOLVER->getSize());
1318 double b = 2.*PI / (right-left) * (symmetric()? 0.5 : 1.0);
1319
1320 assert(TE.rows() == matrixSize());
1321 assert(TH.rows() == matrixSize());
1322
1323 size_t M = TE.cols();
1324 assert(TH.cols() == M);
1325
1326 TempMatrix temp = getTempMatrix();
1327 cmatrix verts(N, M, temp.data());
1328
1329 if (field == FIELD_E) {
1330 if (polarization == E_LONG) {
1331 std::fill_n(verts.data(), verts.rows() * verts.cols(), 0.);
1332 } else if (polarization == E_TRAN) {
1333 if (symmetric()) {
1334 PLASK_OMP_PARALLEL_FOR
1335 for (openmp_size_t m = 0; m < M; m++) {
1336 for (int i = 0; i <= order; ++i) {
1337 dcomplex vert = 0.;
1338 for (int j = 0; j <= order; ++j)
1339 vert += coeff_matrices[l].reyy(i,j) * b*double(j) * TH(j,m);
1340 verts(i,m) = vert / k0;
1341 }
1342 }
1343 } else {
1344 PLASK_OMP_PARALLEL_FOR
1345 for (openmp_size_t m = 0; m < M; m++) {
1346 for (int i = -order; i <= order; ++i) {
1347 dcomplex vert = 0.; // beta is equal to 0
1348 size_t ieh = iEH(i);
1349 for (int j = -order; j <= order; ++j) {
1350 size_t jeh = iEH(j);
1351 vert += coeff_matrices[l].reyy(ieh,jeh) * (b*double(j)-ktran) * TH(jeh,m);
1352 }
1353 verts(ieh,m) = vert / k0;
1354 }
1355 }
1356 }
1357 } else {
1358 if (symmetric()) {
1359 PLASK_OMP_PARALLEL_FOR
1360 for (openmp_size_t m = 0; m < M; m++) {
1361 for (int i = 0; i <= order; ++i) {
1362 dcomplex vert = 0.;
1363 for (int j = 0; j <= order; ++j)
1364 vert -= coeff_matrices[l].reyy(i,j) * (beta * TH(iHx(j),m) + b*double(j) * TH(iHz(j),m));
1365 verts(i,m) = vert / k0;
1366 }
1367 }
1368 } else {
1369 PLASK_OMP_PARALLEL_FOR
1370 for (openmp_size_t m = 0; m < M; m++) {
1371 for (int i = -order; i <= order; ++i) {
1372 dcomplex vert = 0.;
1373 size_t ieh = iEH(i);
1374 for (int j = -order; j <= order; ++j)
1375 vert -= coeff_matrices[l].reyy(ieh,iEH(j)) * (beta * TH(iHx(j),m) + (b*double(j)-ktran) * TH(iHz(j),m));
1376 verts(ieh,m) = vert / k0;
1377 }
1378 }
1379 }
1380 }
1381 } else { // field == FIELD_H
1382 if (polarization == E_TRAN) {
1383 std::fill_n(verts.data(), verts.rows() * verts.cols(), 0.);
1384 } else if (polarization == E_LONG) { // polarization == H_TRAN
1385 if (symmetric()) {
1386 PLASK_OMP_PARALLEL_FOR
1387 for (openmp_size_t m = 0; m < M; m++) {
1388 for (int i = 0; i <= order; ++i) {
1389 dcomplex vert = 0.;
1390 if (periodic)
1391 vert = - (b*double(i)-ktran) * TE(iEH(i),m);
1392 else {
1393 vert = 0.;
1394 for (int j = 0; j <= order; ++j)
1395 vert -= coeff_matrix_rmyy(i,j) * (b*double(j)-ktran) * TE(iEH(j),m);
1396 }
1397 verts(i,m) = vert / k0;
1398 }
1399 }
1400 } else {
1401 PLASK_OMP_PARALLEL_FOR
1402 for (openmp_size_t m = 0; m < M; m++) {
1403 for (int i = -order; i <= order; ++i) {
1404 dcomplex vert = 0.;
1405 size_t ieh = iEH(i);
1406 if (periodic)
1407 vert = - (b*double(i)-ktran) * TE(ieh,m);
1408 else {
1409 vert = 0.;
1410 for (int j = -order; j <= order; ++j) {
1411 size_t jeh = iEH(j);
1412 vert -= coeff_matrix_rmyy(ieh,jeh) * (b*double(j)-ktran) * TE(jeh,m);
1413 }
1414 }
1415 verts(ieh,m) = vert / k0;
1416 }
1417 }
1418 }
1419 } else {
1420 if (symmetric()) {
1421 PLASK_OMP_PARALLEL_FOR
1422 for (openmp_size_t m = 0; m < M; m++) {
1423 for (int i = 0; i <= order; ++i) {
1424 dcomplex vert = 0.;
1425 if (periodic)
1426 vert = (beta * TE(iEx(i),m) - (b*double(i)-ktran) * TE(iEz(i),m));
1427 else {
1428 vert = 0.;
1429 for (int j = 0; j <= order; ++j)
1430 vert += coeff_matrix_rmyy(i,j) * (beta * TE(iEx(j),m) - (b*double(j)-ktran) * TE(iEz(j),m));
1431 }
1432 verts(i,m) = vert / k0;
1433 }
1434 }
1435 } else {
1436 PLASK_OMP_PARALLEL_FOR
1437 for (openmp_size_t m = 0; m < M; m++) {
1438 for (int i = -order; i <= order; ++i) {
1439 dcomplex vert = 0.;
1440 size_t ieh = iEH(i);
1441 if (periodic)
1442 vert = (beta * TE(iEx(i),m) - (b*double(i)-ktran) * TE(iEz(i),m));
1443 else {
1444 vert = 0.;
1445 for (int j = -order; j <= order; ++j)
1446 vert += coeff_matrix_rmyy(ieh,iEH(j)) * (beta * TE(iEx(j),m) - (b*double(j)-ktran) * TE(iEz(j),m));
1447 }
1448 verts(ieh,m) = vert / k0;
1449 }
1450 }
1451 }
1452 }
1453 }
1454
1455 double result = 0.;
1456
1457 if (field == FIELD_E) {
1458 if (symmetric()) {
1459 if (separated()) {
1460 PLASK_OMP_PARALLEL_FOR
1461 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
1462 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
1463 dcomplex sumxz = TE(0,m1) * conj(TE(0,m2)),
1464 sumy = verts(0,m1) * conj(verts(0,m2));
1465 for(size_t i = 1; i != N; ++i) {
1466 sumxz += 2. * TE(i,m1) * conj(TE(i,m2));
1467 sumy += 2. * verts(i,m1) * conj(verts(i,m2));
1468 }
1469 auto vert = vertical(m1, m2);
1470 double res = real(sumxz * vert.first + sumy * vert.second);
1471 if (m2 != m1) res *= 2;
1472 #pragma omp atomic
1473 result += res;
1474 }
1475 }
1476 } else {
1477 PLASK_OMP_PARALLEL_FOR
1478 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
1479 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
1480 size_t ix = iEx(0), iz = iEz(0);
1481 dcomplex sumxz = TE(ix,m1) * conj(TE(ix,m2)) + TE(iz,m1) * conj(TE(iz,m2)),
1482 sumy = verts(0,m1) * conj(verts(0,m2));
1483 for(size_t i = 1; i != N; ++i) {
1484 ix = iEx(i); iz = iEz(i);
1485 sumxz += 2. * (TE(ix,m1) * conj(TE(ix,m2)) + TE(iz,m1) * conj(TE(iz,m2)));
1486 sumy += 2. * verts(i,m1) * conj(verts(i,m2));
1487 }
1488 auto vert = vertical(m1, m2);
1489 double res = real(sumxz * vert.first + sumy * vert.second);
1490 if (m2 != m1) res *= 2;
1491 #pragma omp atomic
1492 result += res;
1493 }
1494 }
1495 }
1496 } else {
1497 if (separated()) {
1498 PLASK_OMP_PARALLEL_FOR
1499 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
1500 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
1501 dcomplex sumxz = 0., sumy = 0.;
1502 for(size_t i = 0; i != N; ++i) {
1503 sumxz += TE(i,m1) * conj(TE(i,m2));
1504 sumy += verts(i,m1) * conj(verts(i,m2));
1505 }
1506 auto vert = vertical(m1, m2);
1507 double res = real(sumxz * vert.first + sumy * vert.second);
1508 if (m2 != m1) res *= 2;
1509 #pragma omp atomic
1510 result += res;
1511 }
1512 }
1513 } else {
1514 PLASK_OMP_PARALLEL_FOR
1515 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
1516 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
1517 dcomplex sumxz = 0., sumy = 0.;
1518 for(size_t i = 0; i != N; ++i) {
1519 size_t ix = iEx(i), iz = iEz(i);
1520 sumxz += TE(ix,m1) * conj(TE(ix,m2)) + TE(iz,m1) * conj(TE(iz,m2));
1521 sumy += verts(i,m1) * conj(verts(i,m2));
1522 }
1523 auto vert = vertical(m1, m2);
1524 double res = real(sumxz * vert.first + sumy * vert.second);
1525 if (m2 != m1) res *= 2;
1526 #pragma omp atomic
1527 result += res;
1528 }
1529 }
1530 }
1531 }
1532 } else { // field == FIELD_H
1533 if (symmetric()) {
1534 if (separated()) {
1535 PLASK_OMP_PARALLEL_FOR
1536 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
1537 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
1538 dcomplex sumxz = TH(0,m1) * conj(TH(0,m2)),
1539 sumy = verts(0,m1) * conj(verts(0,m2));
1540 for(size_t i = 1; i != N; ++i) {
1541 sumxz += 2. * TH(i,m1) * conj(TE(i,m2));
1542 sumy += 2. * verts(i,m1) * conj(verts(i,m2));
1543 }
1544 auto vert = vertical(m1, m2);
1545 double res = real(sumxz * vert.second + sumy * vert.first);
1546 if (m2 != m1) res *= 2;
1547 #pragma omp atomic
1548 result += res;
1549 }
1550 }
1551 } else {
1552 PLASK_OMP_PARALLEL_FOR
1553 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
1554 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
1555 size_t ix = iHx(0), iz = iHz(0);
1556 dcomplex sumxz = TH(ix,m1) * conj(TH(ix,m2)) + TH(iz,m1) * conj(TH(iz,m2)),
1557 sumy = verts(0,m1) * conj(verts(0,m2));
1558 for(size_t i = 1; i != N; ++i) {
1559 ix = iHx(i); iz = iHz(i);
1560 sumxz += 2. * (TH(ix,m1) * conj(TH(ix,m2)) + TE(iz,m1) * conj(TH(iz,m2)));
1561 sumy += 2. * verts(i,m1) * conj(verts(i,m2));
1562 }
1563 auto vert = vertical(m1, m2);
1564 double res = real(sumxz * vert.second + sumy * vert.first);
1565 if (m2 != m1) res *= 2;
1566 #pragma omp atomic
1567 result += res;
1568 }
1569 }
1570 }
1571 } else {
1572 if (separated()) {
1573 PLASK_OMP_PARALLEL_FOR
1574 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
1575 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
1576 dcomplex sumxz = 0., sumy = 0.;
1577 for(size_t i = 0; i != N; ++i) {
1578 sumxz += TH(i,m1) * conj(TH(i,m2));
1579 sumy += verts(i,m1) * conj(verts(i,m2));
1580 }
1581 auto vert = vertical(m1, m2);
1582 double res = real(sumxz * vert.second + sumy * vert.first);
1583 if (m2 != m1) res *= 2;
1584 #pragma omp atomic
1585 result += res;
1586 }
1587 }
1588 } else {
1589 PLASK_OMP_PARALLEL_FOR
1590 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
1591 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
1592 dcomplex sumxz = 0., sumy = 0.;
1593 for(size_t i = 0; i != N; ++i) {
1594 size_t ix = iHx(i), iz = iHz(i);
1595 sumxz += TH(ix,m1) * conj(TH(ix,m2)) + TH(iz,m1) * conj(TH(iz,m2));
1596 sumy += verts(i,m1) * conj(verts(i,m2));
1597 }
1598 auto vert = vertical(m1, m2);
1599 double res = real(sumxz * vert.second + sumy * vert.first);
1600 if (m2 != m1) res *= 2;
1601 #pragma omp atomic
1602 result += res;
1603 }
1604 }
1605 }
1606 }
1607 }
1608
1609 const double L = (right-left) * (symmetric()? 2. : 1.);
1610 return 0.5 * result * L;
1611}
1612
1613
1614}}} // namespace plask