PLaSK library

Loading...
Searching...
No Matches
expansioncyl.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 "expansioncyl.hpp"
15#include "solvercyl.hpp"
16#include "zeros-data.hpp"
17
18#include "expansioncyl-fini.hpp"
19#include "expansioncyl-infini.hpp"
20
21#include "../gauss_legendre.hpp"
22
23#include "besselj.hpp"
24
25#include <boost/math/special_functions/legendre.hpp>
26 using boost::math::legendre_p;
27
28#define SOLVER static_cast<BesselSolverCyl*>(solver)
29
30 namespace plask { namespace optical { namespace modal {
31
32 ExpansionBessel::ExpansionBessel(BesselSolverCyl* solver) : Expansion(solver), m(1), initialized(false), m_changed(true) {}
33
34 size_t ExpansionBessel::matrixSize() const { return 2 * SOLVER->size; }
35
36 void ExpansionBessel::init1() {
37 // Initialize segments
38 if (SOLVER->mesh)
39 rbounds = OrderedAxis(*SOLVER->getMesh());
40 else
41 rbounds = std::move(*makeGeometryGrid1D(SOLVER->getGeometry()));
42 OrderedAxis::WarningOff nowarn_rbounds(rbounds);
43 rbounds.addPoint(0.);
44 size_t nseg = rbounds.size() - 1;
45 if (dynamic_cast<ExpansionBesselFini*>(this)) {
46 if (SOLVER->pml.dist > 0.) rbounds.addPoint(rbounds[nseg++] + SOLVER->pml.dist);
47 if (SOLVER->pml.size > 0.) rbounds.addPoint(rbounds[nseg++] + SOLVER->pml.size);
48 }
49 segments.resize(nseg);
50 double a, b = 0.;
51 for (size_t i = 0; i < nseg; ++i) {
52 a = b;
53 b = rbounds[i + 1];
54 segments[i].Z = 0.5 * (a + b);
55 segments[i].D = 0.5 * (b - a);
56 }
57
58 diagonals.assign(solver->lcount, false);
59 initialized = true;
60 m_changed = true;
61}
62
63 void ExpansionBessel::reset() {
64 layers_integrals.clear();
65 segments.clear();
66 kpts.clear();
67 initialized = false;
68 mesh.reset();
69 temporary.reset();
70}
71
72 void ExpansionBessel::init3() {
73 size_t nseg = rbounds.size() - 1;
74
75 // Estimate necessary number of integration points
76 double k = kpts[kpts.size() - 1];
77
78 double expected = rbounds[rbounds.size() - 1] * cyl_bessel_j(m + 1, k);
79 expected = 0.5 * expected * expected;
80
81 k /= rbounds[rbounds.size() - 1];
82
83 double max_error = SOLVER->integral_error * expected / double(nseg);
84 double error = 0.;
85
86 std::deque<std::vector<double>> abscissae_cache;
87 std::deque<DataVector<double>> weights_cache;
88
89 auto raxis = plask::make_shared<OrderedAxis>();
90 OrderedAxis::WarningOff nowarn_raxis(raxis);
91
92 double expcts = 0.;
93 for (size_t i = 0; i < nseg; ++i) {
94 double b = rbounds[i + 1];
95
96 // expected value is the second Lommel's integral
97 double expct = expcts;
98 expcts = cyl_bessel_j(m, k * b);
99 expcts = 0.5 * b * b * (expcts * expcts - cyl_bessel_j(m - 1, k * b) * cyl_bessel_j(m + 1, k * b));
100 expct = expcts - expct;
101
102 double err = 2 * max_error;
103 std::vector<double> points;
104 size_t j, n = 0;
105 double sum;
106 for (j = 0; err > max_error && n <= SOLVER->max_integration_points; ++j) {
107 n = 4 * (j + 1) - 1;
108 if (j == abscissae_cache.size()) {
109 abscissae_cache.push_back(std::vector<double>());
110 weights_cache.push_back(DataVector<double>());
111 gaussLegendre(n, abscissae_cache.back(), weights_cache.back());
112 }
113 assert(j < abscissae_cache.size());
114 assert(j < weights_cache.size());
115 const std::vector<double>& abscissae = abscissae_cache[j];
116 points.clear();
117 points.reserve(abscissae.size());
118 sum = 0.;
119 for (size_t a = 0; a != abscissae.size(); ++a) {
120 double r = segments[i].Z + segments[i].D * abscissae[a];
121 double Jm = cyl_bessel_j(m, k * r);
122 sum += weights_cache[j][a] * r * Jm * Jm;
123 points.push_back(r);
124 }
125 sum *= segments[i].D;
126 err = abs(sum - expct);
127 }
128 error += err;
129 raxis->addOrderedPoints(points.begin(), points.end());
130 segments[i].weights = weights_cache[j - 1];
131 }
132
133 SOLVER->writelog(LOG_DETAIL, "Sampling structure in {:d} points (error: {:g}/{:g})", raxis->size(), error / expected,
134 SOLVER->integral_error);
135
136 // Allocate memory for integrals
137 size_t nlayers = solver->lcount;
138 layers_integrals.resize(nlayers);
139
140 mesh = plask::make_shared<RectangularMesh<2>>(raxis, solver->verts, RectangularMesh<2>::ORDER_01);
141
142 m_changed = false;
143}
144
145 void ExpansionBessel::beforeLayersIntegrals(dcomplex lam, dcomplex glam) {
146 if (m_changed) init2();
147 SOLVER->prepareExpansionIntegrals(this, mesh, lam, glam);
148}
149
150 void ExpansionBessel::beforeGetEpsilon() {
151 double lambda = real(2e3 * PI / k0);
152 double lam, glam;
153 if (!isnan(lam0)) {
154 lam = lam0;
155 glam = (solver->always_recompute_gain) ? lambda : lam;
156 } else {
157 lam = glam = lambda;
158 }
159 beforeLayersIntegrals(lam, glam);
160}
161
162 void ExpansionBessel::afterGetEpsilon() { afterLayersIntegrals(); }
163
164 Tensor3<dcomplex> ExpansionBessel::getEps(size_t layer, size_t ri, double r, double matz, double lam, double glam) {
165 Tensor3<dcomplex> eps;
166 std::set<std::string> roles;
167 if (epsilon_connected && solver->lcomputed[layer] || gain_connected && solver->lgained[layer])
168 roles = SOLVER->getGeometry()->getRolesAt(vec(r, matz));
169 bool computed = epsilon_connected && solver->lcomputed[layer] && roles.find("inEpsilon") != roles.end();
170 if (computed) {
171 eps = Zero<Tensor3<dcomplex>>();
172 double W = 0.;
173 for (size_t k = 0, v = ri * solver->verts->size(); k != mesh->vert()->size(); ++v, ++k) {
174 if (solver->stack[k] == layer) {
175 if (isnan(epsilons[v]))
176 throw BadInput(solver->getId(), "complex permittivity tensor got from inEpsilon is NaN at {}", mesh->at(v));
177 double w = (k == 0 || k == mesh->vert()->size() - 1) ? 1e-6 : solver->vbounds->at(k) - solver->vbounds->at(k - 1);
178 eps += w * epsilons[v];
179 W += w;
180 }
181 }
182 eps /= W;
183 } else {
184 OmpLockGuard lock; // this must be declared before `material` to guard its destruction
185 auto material = SOLVER->getGeometry()->getMaterial(vec(r, matz));
186 lock = material->lock();
187 double T, cc;
188 std::tie(T, cc) = getTC(layer, ri);
189 eps = material->Eps(lam, T, cc);
190 if (isnan(eps))
191 throw BadInput(solver->getId(), "complex permittivity tensor (Eps) for {} is NaN at lam={}nm, T={}K, n={}/cm3",
192 material->name(), lam, T, cc);
193 }
194 if (!is_zero(eps.c00 - eps.c11) || eps.c01 != 0. || eps.c02 != 0. || eps.c10 != 0. || eps.c12 != 0. || eps.c20 != 0. ||
195 eps.c21 != 0.)
196 throw BadInput(solver->getId(), "lateral anisotropy not allowed for this solver");
197 if (!computed && gain_connected && solver->lgained[layer]) {
198 if (roles.find("QW") != roles.end() || roles.find("QD") != roles.end() || roles.find("gain") != roles.end()) {
199 Tensor2<double> g = 0.;
200 double W = 0.;
201 for (size_t k = 0, v = ri * solver->verts->size(); k != mesh->vert()->size(); ++v, ++k) {
202 if (solver->stack[k] == layer) {
203 double w =
204 (k == 0 || k == mesh->vert()->size() - 1) ? 1e-6 : solver->vbounds->at(k) - solver->vbounds->at(k - 1);
205 g += w * gain[v];
206 W += w;
207 }
208 }
209 Tensor2<double> ni = glam * g / W * (0.25e-7 / PI);
210 double n00 = sqrt(eps.c00).real(), n22 = sqrt(eps.c22).real();
211 eps.c00 = eps.c11 = dcomplex(n00 * n00 - ni.c00 * ni.c00, 2 * n00 * ni.c00);
212 eps.c22 = dcomplex(n22 * n22 - ni.c11 * ni.c11, 2 * n22 * ni.c11);
213 }
214 }
215 return eps;
216}
217
218 void ExpansionBessel::layerIntegrals(size_t layer, double lam, double glam) {
219 if (isnan(real(k0)) || isnan(imag(k0))) throw BadInput(SOLVER->getId(), "no wavelength specified");
220
221 auto geometry = SOLVER->getGeometry();
222
223 auto raxis = mesh->tran();
224
225#if defined(OPENMP_FOUND) // && !defined(NDEBUG)
226 SOLVER->writelog(LOG_DETAIL, "Computing integrals for layer {:d}/{:d} with {} rule in thread {:d}", layer, solver->lcount,
227 SOLVER->ruleName(), omp_get_thread_num());
228#else
229 SOLVER->writelog(LOG_DETAIL, "Computing integrals for layer {:d}/{:d} with {} rule", layer, solver->lcount, SOLVER->ruleName());
230#endif
231
232 if (isnan(lam)) throw BadInput(SOLVER->getId(), "no wavelength given: specify 'lam' or 'lam0'");
233
234 size_t nr = raxis->size(), N = SOLVER->size;
235
236 if (epsilon_connected && solver->lcomputed[layer]) {
237 SOLVER->writelog(LOG_DEBUG, "Layer {:d} takes some materials parameters from inEpsilon", layer);
238 if (isnan(glam)) glam = lam;
239 }
240 if (gain_connected && solver->lgained[layer]) {
241 SOLVER->writelog(LOG_DEBUG, "Layer {:d} has gain", layer);
242 if (isnan(glam)) glam = lam;
243 }
244
245 double matz;
246 for (size_t i = 0; i != solver->stack.size(); ++i) {
247 if (solver->stack[i] == layer) {
248 matz = solver->verts->at(i);
249 break;
250 }
251 }
252
253 // For checking if the layer is uniform
254 diagonals[layer] = true;
255
256 bool finite = !dynamic_cast<ExpansionBesselInfini*>(this);
257
258 size_t pmli = raxis->size();
259 double pmlr;
260 dcomplex epsp0, epsr0, epsz0;
261 if (finite) {
262 if (SOLVER->pml.size > 0. && SOLVER->pml.factor != 1.) {
263 size_t pmlseg = segments.size() - 1;
264 pmli -= segments[pmlseg].weights.size();
265 pmlr = rbounds[pmlseg];
266 }
267 } else {
268 Tensor3<dcomplex> eps0 = getEps(layer, mesh->tran()->size() - 1, rbounds[rbounds.size() - 1] + 0.001, matz, lam, glam);
269 epsp0 = eps0.c00;
270 if (SOLVER->rule != BesselSolverCyl::RULE_OLD) {
271 epsr0 = (SOLVER->rule != BesselSolverCyl::RULE_DIRECT) ? 1. / eps0.c11 : eps0.c11;
272 epsz0 = eps0.c22;
273 } else {
274 epsr0 = eps0.c11;
275 epsz0 = 1. / eps0.c22;
276 }
277 if (abs(epsp0.imag()) < SMALL) epsp0.imag(0.);
278 if (abs(epsr0.imag()) < SMALL) epsr0.imag(0.);
279 if (abs(epsz0.imag()) < SMALL) epsz0.imag(0.);
280 writelog(LOG_DEBUG, "Reference refractive index for layer {} is {} / {}", layer, str(sqrt(epsp0)),
281 str(sqrt((SOLVER->rule != BesselSolverCyl::RULE_OLD) ? epsz0 : (1. / epsz0))));
282 }
283
284 aligned_unique_ptr<dcomplex> epsp_data(aligned_malloc<dcomplex>(nr));
285 aligned_unique_ptr<dcomplex> epsr_data(aligned_malloc<dcomplex>(nr));
286 aligned_unique_ptr<dcomplex> epsz_data(aligned_malloc<dcomplex>(nr));
287
288 // Compute integrals
289 for (size_t ri = 0, wi = 0, seg = 0, nw = segments[0].weights.size(); ri != nr; ++ri, ++wi) {
290 if (wi == nw) {
291 nw = segments[++seg].weights.size();
292 wi = 0;
293 }
294 double r = raxis->at(ri);
295 double w = segments[seg].weights[wi] * segments[seg].D;
296
297 Tensor3<dcomplex> eps = getEps(layer, ri, r, matz, lam, glam);
298 if (ri >= pmli) {
299 dcomplex f = 1. + (SOLVER->pml.factor - 1.) * pow((r - pmlr) / SOLVER->pml.size, SOLVER->pml.order);
300 eps.c00 *= f;
301 eps.c11 /= f;
302 eps.c22 *= f;
303 }
304 dcomplex epsp = eps.c00, epsr, epsz;
305 if (SOLVER->rule != BesselSolverCyl::RULE_OLD) {
306 epsr = (SOLVER->rule != BesselSolverCyl::RULE_DIRECT) ? 1. / eps.c11 : eps.c11;
307 epsz = eps.c22;
308 } else {
309 epsr = eps.c11;
310 epsz = 1. / eps.c22;
311 }
312
313 if (finite) {
314 if (ri == 0) {
315 epsp0 = epsp;
316 epsr0 = epsr;
317 epsz0 = epsz;
318 } else {
319 if (!is_zero(epsp - epsp0) || !is_zero(epsr - epsr0) || !is_zero(epsz - epsz0)) diagonals[layer] = false;
320 }
321 } else {
322 epsp -= epsp0;
323 epsr -= epsr0;
324 epsz -= epsz0;
325 if (!is_zero(epsp) || !is_zero(epsr) || !is_zero(epsz)) diagonals[layer] = false;
326 }
327
328 epsp_data.get()[ri] = epsp * w;
329 epsr_data.get()[ri] = epsr * w;
330 epsz_data.get()[ri] = epsz * w;
331 }
332
333 if (diagonals[layer]) {
334 Integrals& integrals = layers_integrals[layer];
335 SOLVER->writelog(LOG_DETAIL, "Layer {0} is uniform", layer);
336 integrals.reset(N);
337 zero_matrix(integrals.V_k);
338 zero_matrix(integrals.Tss);
339 zero_matrix(integrals.Tsp);
340 zero_matrix(integrals.Tps);
341 zero_matrix(integrals.Tpp);
342 double ib = 1. / rbounds[rbounds.size() - 1];
343 for (size_t i = 0; i < N; ++i) {
344 double k = kpts[i] * ib;
345 if (SOLVER->rule != BesselSolverCyl::RULE_OLD)
346 integrals.V_k(i, i) = k / epsz0;
347 else
348 integrals.V_k(i, i) = k * epsz0;
349 // integrals.Tss(i,i) = integrals.Tpp(i,i) = 1. / iepsr0 + epsp0;
350 // integrals.Tsp(i,i) = integrals.Tps(i,i) = 1. / iepsr0 - epsp0;
351 integrals.Tss(i, i) = integrals.Tpp(i, i) = 2. * epsp0;
352 }
353 } else {
354 integrateParams(layers_integrals[layer], epsp_data.get(), epsr_data.get(), epsz_data.get(), epsp0, epsr0, epsz0);
355 }
356}
357
358#ifndef NDEBUG
359 cmatrix ExpansionBessel::epsV_k(size_t layer) {
360 size_t N = SOLVER->size;
361 cmatrix result(N, N, 0.);
362 for (size_t i = 0; i != N; ++i)
363 for (size_t j = 0; j != N; ++j) result(i, j) = layers_integrals[layer].V_k(i, j);
364 return result;
365}
366 cmatrix ExpansionBessel::epsTss(size_t layer) {
367 size_t N = SOLVER->size;
368 cmatrix result(N, N, 0.);
369 for (size_t i = 0; i != N; ++i)
370 for (size_t j = 0; j != N; ++j) result(i, j) = layers_integrals[layer].Tss(i, j);
371 return result;
372}
373 cmatrix ExpansionBessel::epsTsp(size_t layer) {
374 size_t N = SOLVER->size;
375 cmatrix result(N, N, 0.);
376 for (size_t i = 0; i != N; ++i)
377 for (size_t j = 0; j != N; ++j) result(i, j) = layers_integrals[layer].Tsp(i, j);
378 return result;
379}
380 cmatrix ExpansionBessel::epsTps(size_t layer) {
381 size_t N = SOLVER->size;
382 cmatrix result(N, N, 0.);
383 for (size_t i = 0; i != N; ++i)
384 for (size_t j = 0; j != N; ++j) result(i, j) = layers_integrals[layer].Tps(i, j);
385 return result;
386}
387 cmatrix ExpansionBessel::epsTpp(size_t layer) {
388 size_t N = SOLVER->size;
389 cmatrix result(N, N, 0.);
390 for (size_t i = 0; i != N; ++i)
391 for (size_t j = 0; j != N; ++j) result(i, j) = layers_integrals[layer].Tpp(i, j);
392 return result;
393}
394#endif
395
396 void ExpansionBessel::prepareField() {
397 if (field_interpolation == INTERPOLATION_DEFAULT) field_interpolation = INTERPOLATION_NEAREST;
398}
399
400 void ExpansionBessel::cleanupField() {}
401
402 LazyData<Vec<3, dcomplex>> ExpansionBessel::getField(size_t layer,
403 const shared_ptr<const typename LevelsAdapter::Level>& level,
404 const cvector& E,
405 const cvector& H) {
406 size_t N = SOLVER->size;
407
408 assert(dynamic_pointer_cast<const MeshD<2>>(level->mesh()));
409 auto dest_mesh = static_pointer_cast<const MeshD<2>>(level->mesh());
410 double ib = 1. / rbounds[rbounds.size() - 1];
411 const dcomplex fz = -I / k0;
412
413 auto src_mesh =
414 plask::make_shared<RectangularMesh<2>>(mesh->tran(), plask::make_shared<RegularAxis>(level->vpos(), level->vpos(), 1));
415
416 if (which_field == FIELD_E) {
417 cvector Ez(N);
418 {
419 cvector Dz(N);
420 for (size_t j = 0; j != N; ++j) {
421 size_t js = idxs(j), jp = idxp(j);
422 Dz[j] = fz * (H[js] + H[jp]);
423 }
424 mult_matrix_by_vector(layers_integrals[layer].V_k, Dz, Ez);
425 }
426 return LazyData<Vec<3, dcomplex>>(dest_mesh->size(), [=](size_t i) -> Vec<3, dcomplex> {
427 double r = dest_mesh->at(i)[0];
428 Vec<3, dcomplex> result{0., 0., 0.};
429 for (size_t j = 0; j != N; ++j) {
430 double k = kpts[j] * ib;
431 double kr = k * r;
432 double Jm = cyl_bessel_j(m - 1, kr), Jp = cyl_bessel_j(m + 1, kr), J = cyl_bessel_j(m, kr);
433 size_t js = idxs(j), jp = idxp(j);
434 result.c0 += Jm * E[js] - Jp * E[jp]; // E_p
435 result.c1 += Jm * E[js] + Jp * E[jp]; // E_r
436 result.c2 += J * Ez[j]; // E_z
437 }
438 return result;
439 });
440 } else { // which_field == FIELD_H
441 double r0 = (SOLVER->pml.size > 0. && SOLVER->pml.factor != 1.) ? rbounds[rbounds.size() - 1] : INFINITY;
442 return LazyData<Vec<3, dcomplex>>(dest_mesh->size(), [=](size_t i) -> Vec<3, dcomplex> {
443 double r = dest_mesh->at(i)[0];
444 dcomplex imu = 1.;
445 if (r > r0) imu = 1. / (1. + (SOLVER->pml.factor - 1.) * pow((r - r0) / SOLVER->pml.size, SOLVER->pml.order));
446 Vec<3, dcomplex> result{0., 0., 0.};
447 for (size_t j = 0; j != N; ++j) {
448 double k = kpts[j] * ib;
449 double kr = k * r;
450 double Jm = cyl_bessel_j(m - 1, kr), Jp = cyl_bessel_j(m + 1, kr), J = cyl_bessel_j(m, kr);
451 size_t js = idxs(j), jp = idxp(j);
452 result.c0 -= Jm * H[js] - Jp * H[jp]; // H_p
453 result.c1 += Jm * H[js] + Jp * H[jp]; // H_r
454 result.c2 += fz * k * imu * J * (E[js] + E[jp]); // H_z
455 }
456 return result;
457 });
458 }
459}
460
461 LazyData<Tensor3<dcomplex>> ExpansionBessel::getMaterialEps(size_t layer,
462 const shared_ptr<const typename LevelsAdapter::Level>& level,
463 InterpolationMethod interp) {
464 if (interp == INTERPOLATION_DEFAULT) interp = INTERPOLATION_NEAREST;
465
466 assert(dynamic_pointer_cast<const MeshD<2>>(level->mesh()));
467 auto dest_mesh = static_pointer_cast<const MeshD<2>>(level->mesh());
468
469 double lam, glam;
470 if (!isnan(lam0)) {
471 lam = lam0;
472 glam = (solver->always_recompute_gain) ? real(2e3 * PI / k0) : lam;
473 } else {
474 lam = glam = real(2e3 * PI / k0);
475 }
476
477 auto raxis = mesh->tran();
478
479 DataVector<Tensor3<dcomplex>> eps(raxis->size());
480 for (size_t i = 0; i != eps.size(); ++i) {
481 Tensor3<dcomplex> eps = getEps(layer, i, raxis->at(i), level->vpos(), lam, glam);
482 }
483
484 auto src_mesh =
485 plask::make_shared<RectangularMesh<2>>(mesh->tran(), plask::make_shared<RegularAxis>(level->vpos(), level->vpos(), 1));
486 return interpolate(
487 src_mesh, eps, dest_mesh, interp,
488 InterpolationFlags(SOLVER->getGeometry(), InterpolationFlags::Symmetry::POSITIVE, InterpolationFlags::Symmetry::NO));
489}
490
491 double ExpansionBessel::integratePoyntingVert(const cvector& E, const cvector& H) {
492 double result = 0.;
493 for (size_t i = 0, N = SOLVER->size; i < N; ++i) {
494 size_t is = idxs(i);
495 size_t ip = idxp(i);
496 result += real(-E[is] * conj(H[is]) + E[ip] * conj(H[ip])) * fieldFactor(i);
497 }
498 return 4e-12 * PI * result; // µm² -> m²
499}
500
501 double ExpansionBessel::integrateField(WhichField field,
502 size_t layer,
503 const cmatrix& TE,
504 const cmatrix& TH,
505 const std::function<std::pair<dcomplex, dcomplex>(size_t, size_t)>& vertical) {
506 assert(TE.rows() == matrixSize());
507 assert(TH.rows() == matrixSize());
508
509 size_t M = TE.cols();
510 assert(TH.cols() == M);
511
512 size_t N = SOLVER->size;
513
514 TempMatrix temp = getTempMatrix();
515 cmatrix Fz(N, M, temp.data()), DBz(N, M, temp.data() + N * M);
516
517 double R = rbounds[rbounds.size() - 1];
518 double fz = 0.5 / real(k0 * conj(k0));
519
520 if (field == FIELD_E) {
521 PLASK_OMP_PARALLEL_FOR
522 for (openmp_size_t m = 0; m < M; m++) {
523 cvector Ez(N), Dz(N);
524 for (size_t j = 0; j != N; ++j) {
525 size_t js = idxs(j), jp = idxp(j);
526 DBz(j, m) = TH(js, m) + TH(jp, m);
527 }
528 }
529 mult_matrix_by_matrix(layers_integrals[layer].V_k, DBz, Fz);
530 } else {
531 PLASK_OMP_PARALLEL_FOR
532 for (openmp_size_t m = 0; m < M; m++) {
533 for (size_t j = 0; j != N; ++j) {
534 size_t js = idxs(j), jp = idxp(j);
535 DBz(j, m) = TE(js, m) + TE(jp, m);
536 }
537 }
538 Fz = getHzMatrix(DBz, Fz);
539 }
540
541 double result = 0.;
542
543 if (field == FIELD_E) {
544 PLASK_OMP_PARALLEL_FOR
545 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
546 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
547 dcomplex resxy = 0., resz = 0.;
548 for (size_t i = 0, N = SOLVER->size; i < N; ++i) {
549 double eta = fieldFactor(i);
550 size_t is = idxs(i);
551 size_t ip = idxp(i);
552 resxy += (TE(is, m1) * conj(TE(is, m2)) + TE(ip, m1) * conj(TE(ip, m2))) * eta;
553 resz += Fz(i, m1) * conj(Fz(i, m2)) * eta;
554 }
555 if (!(is_zero(resxy) && is_zero(resz))) {
556 auto vert = vertical(m1, m2);
557 double res = real(resxy * vert.first + fz * resz * vert.second);
558 if (m2 != m1) res *= 2;
559#pragma omp atomic
560 result += res;
561 }
562 }
563 }
564 } else {
565 PLASK_OMP_PARALLEL_FOR
566 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
567 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
568 dcomplex resxy = 0., resz = 0.;
569 for (size_t i = 0, N = SOLVER->size; i < N; ++i) {
570 double eta = fieldFactor(i);
571 size_t is = idxs(i);
572 size_t ip = idxp(i);
573 resxy += (TH(is, m1) * conj(TH(is, m2)) + TH(ip, m1) * conj(TH(ip, m2))) * eta;
574 resz += Fz(i, m1) * conj(Fz(i, m2)) * eta;
575 }
576 if (!(is_zero(resxy) && is_zero(resz))) {
577 auto vert = vertical(m1, m2);
578 double res = real(resxy * vert.second + fz * resz * vert.first);
579 if (m2 != m1) res *= 2;
580#pragma omp atomic
581 result += res;
582 }
583 }
584 }
585 }
586 return 2 * PI * result;
587}
588
589}}} // namespace plask::optical::modal