PLaSK library

Loading...
Searching...
No Matches
efm.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 "efm.hpp"
15#include "gauss_matrix.hpp"
16#include "patterson.hpp"
17
18 using plask::dcomplex;
19
20 namespace plask { namespace optical { namespace effective {
21
22#define DLAM 1e-3
23
24 EffectiveFrequencyCyl::EffectiveFrequencyCyl(const std::string& name)
25 : SolverWithMesh<Geometry2DCylindrical, RectangularMesh<2>>(name),
26 log_value(dataLog<dcomplex, dcomplex>("radial", "lam", "det")),
27 emission(TOP),
28 rstripe(-1),
29 determinant(DETERMINANT_OUTWARDS),
30 perr(1e-3),
31 k0(NAN),
32 vlam(0.),
33 outWavelength(this, &EffectiveFrequencyCyl::getWavelength, &EffectiveFrequencyCyl::nmodes),
34 outLoss(this, &EffectiveFrequencyCyl::getModalLoss, &EffectiveFrequencyCyl::nmodes),
35 outLightMagnitude(this, &EffectiveFrequencyCyl::getLightMagnitude, &EffectiveFrequencyCyl::nmodes),
36 outLightE(this, &EffectiveFrequencyCyl::getElectricField, &EffectiveFrequencyCyl::nmodes),
37 outRefractiveIndex(this, &EffectiveFrequencyCyl::getRefractiveIndex),
38 outHeat(this, &EffectiveFrequencyCyl::getHeat),
39 asymptotic(false) {
40 inTemperature = 300.;
41 root.tolx = 1.0e-6;
42 root.tolf_min = 1.0e-7;
43 root.tolf_max = 2.0e-5;
44 root.maxiter = 500;
45 root.method = RootDigger::ROOT_MULLER;
46 stripe_root.tolx = 1.0e-6;
47 stripe_root.tolf_min = 1.0e-7;
48 stripe_root.tolf_max = 1.0e-5;
49 stripe_root.maxiter = 500;
50 stripe_root.method = RootDigger::ROOT_MULLER;
51 inTemperature.changedConnectMethod(this, &EffectiveFrequencyCyl::onInputChange);
52 inGain.changedConnectMethod(this, &EffectiveFrequencyCyl::onInputChange);
53 inCarriersConcentration.changedConnectMethod(this, &EffectiveFrequencyCyl::onInputChange);
54}
55
56 void EffectiveFrequencyCyl::loadConfiguration(XMLReader& reader, Manager& manager) {
57 while (reader.requireTagOrEnd()) {
58 std::string param = reader.getNodeName();
59 if (param == "mode") {
60 // m = reader.getAttribute<unsigned short>("m", m);
61 auto alam0 = reader.getAttribute<double>("lam0");
62 auto ak0 = reader.getAttribute<double>("k0");
63 if (alam0) {
64 if (ak0) throw XMLConflictingAttributesException(reader, "k0", "lam0");
65 k0 = 2e3 * PI / *alam0;
66 } else if (ak0)
67 k0 = *ak0;
68 emission = reader.enumAttribute<Emission>("emission").value("top", TOP).value("bottom", BOTTOM).get(emission);
69 vlam = reader.getAttribute<double>("vlam", real(vlam));
70 if (reader.hasAttribute("vat")) {
71 std::string str = reader.requireAttribute("vat");
72 if (str == "all" || str == "none")
73 rstripe = -1;
74 else
75 try {
76 setStripeR(boost::lexical_cast<double>(str));
77 } catch (boost::bad_lexical_cast&) {
78 throw XMLBadAttrException(reader, "vat", str);
79 }
80 }
81 asymptotic = reader.getAttribute<bool>("asymptotic", asymptotic);
82 reader.requireTagEnd();
83 } else if (param == "root") {
84 determinant = reader.enumAttribute<Determinant>("determinant")
85 .value("full", DETERMINANT_FULL)
86 .value("outwards", DETERMINANT_OUTWARDS)
87 .value("inwards", DETERMINANT_INWARDS)
88 .value("transfer", DETERMINANT_INWARDS)
89 .get(determinant);
90 RootDigger::readRootDiggerConfig(reader, root);
91 } else if (param == "stripe-root") {
92 RootDigger::readRootDiggerConfig(reader, stripe_root);
93 } else if (param == "mesh") {
94 auto name = reader.getAttribute("ref");
95 if (!name)
96 name.reset(reader.requireTextInCurrentTag());
97 else
98 reader.requireTagEnd();
99 auto found = manager.meshes.find(*name);
100 if (found != manager.meshes.end()) {
101 auto mesh1 = dynamic_pointer_cast<MeshAxis>(found->second);
102 auto mesh2 = dynamic_pointer_cast<RectangularMesh<2>>(found->second);
103 if (mesh1)
104 this->setHorizontalMesh(mesh1);
105 else if (mesh2)
106 this->setMesh(mesh2);
107 else {
108 auto generator1 = dynamic_pointer_cast<MeshGeneratorD<1>>(found->second);
109 auto generator2 = dynamic_pointer_cast<MeshGeneratorD<2>>(found->second);
110 if (generator1)
111 this->setMesh(plask::make_shared<RectangularMesh2DFrom1DGenerator>(generator1));
112 else if (generator2)
113 this->setMesh(generator2);
114 else
115 throw BadInput(this->getId(), "mesh or generator '{0}' of wrong type", *name);
116 }
117 }
118 } else
119 parseStandardConfiguration(reader, manager, "<geometry>, <mesh>, <mode>, <root>, <stripe-root>, or <outer>");
120 }
121}
122
123 size_t EffectiveFrequencyCyl::findMode(dcomplex lambda, int m) {
124 writelog(LOG_INFO, "Searching for the mode starting from wavelength = {0}", str(lambda));
125 if (isnan(k0.real())) throw BadInput(getId(), "no reference wavelength `lam0` specified");
126 stageOne();
127 Mode mode(this, m);
128 mode.lam = RootDigger::get(
129 this, [this, &mode](const dcomplex& lam) { return this->detS(lam, mode); }, log_value, root)
130 ->find(lambda);
131 return insertMode(mode);
132}
133
134std::vector<size_t> EffectiveFrequencyCyl::findModes(dcomplex lambda1,
135 dcomplex lambda2,
136 int m,
137 size_t resteps,
138 size_t imsteps,
139 dcomplex eps) {
140 if (isnan(k0.real())) throw BadInput(getId(), "no reference wavelength `lam0` specified");
141 if (isnan(k0.real())) throw BadInput(getId(), "no reference wavelength `lam0` specified");
142 stageOne();
143
144 if ((real(lambda1) == 0. && real(lambda2) != 0.) || (real(lambda1) != 0. && real(lambda2) == 0.))
145 throw BadInput(getId(), "bad area to browse specified");
146
147 dcomplex lam0 = lambda1;
148 dcomplex lam1 = lambda2;
149
150 if (eps.imag() == 0.) eps.imag(eps.real());
151
152 if (real(eps) <= 0. || imag(eps) <= 0.) throw BadInput(this->getId(), "bad precision specified");
153
154 double re0 = real(lam0), im0 = imag(lam0);
155 double re1 = real(lam1), im1 = imag(lam1);
156 if (re0 > re1) std::swap(re0, re1);
157 if (im0 > im1) std::swap(im0, im1);
158
159 if (real(lambda1) == 0. && real(lambda2) == 0.) {
160 re0 = 1e30;
161 re1 = -1e30;
162 for (size_t i = 0; i != rsize; ++i) {
163 dcomplex v = veffs[i];
164 if (v.real() < re0) re0 = v.real();
165 if (v.real() > re1) re1 = v.real();
166 }
167 }
168 if (imag(lambda1) == 0. && imag(lambda2) == 0.) {
169 im0 = 1e30;
170 im1 = -1e30;
171 for (size_t i = 0; i != rsize; ++i) {
172 dcomplex v = veffs[i];
173 if (v.imag() < im0) im0 = v.imag();
174 if (v.imag() > im1) im1 = v.imag();
175 }
176 }
177 lam0 = 1.000001 * dcomplex(re0, im0);
178 lam1 = 0.999999 * dcomplex(re1, im1);
179
180 Mode mode(this, m);
181 auto results = findZeros(
182 this, [this, &mode](dcomplex lam) { return this->detS(lam, mode); }, lam0, lam1, resteps, imsteps, eps);
183
184 std::vector<size_t> idx(results.size());
185
186 if (results.size() != 0) {
187 log_value.resetCounter();
188 auto refine = RootDigger::get(
189 this, [this, &mode](const dcomplex& lam) { return this->detS(lam, mode); }, log_value, root);
190 std::string msg = "Found modes at: ";
191 for (auto& zz : results) {
192 dcomplex z;
193 try {
194 z = refine->find(0.5 * (zz.first + zz.second));
195 } catch (ComputationError&) {
196 continue;
197 }
198 mode.lam = z;
199 idx.push_back(insertMode(mode));
200 msg += str(z) + ", ";
201 }
202 writelog(LOG_RESULT, msg.substr(0, msg.length() - 2));
203 } else
204 writelog(LOG_RESULT, "Did not find any modes");
205
206 return idx;
207}
208
209 size_t EffectiveFrequencyCyl::setMode(dcomplex clambda, int m) {
210 if (isnan(k0.real())) throw BadInput(getId(), "no reference wavelength `lam0` specified");
211 stageOne();
212 Mode mode(this, m);
213 mode.lam = clambda;
214 double det = abs(detS(mode.lam, mode));
215 if (det > root.tolf_max) writelog(LOG_WARNING, "Provided wavelength does not correspond to any mode (det = {0})", det);
216 writelog(LOG_INFO, "Setting mode at {0}", str(clambda));
217 return insertMode(mode);
218}
219
220 void EffectiveFrequencyCyl::onInitialize() {
221 if (!geometry) throw NoGeometryException(getId());
222
223 // Set default mesh
224 if (!mesh) setSimpleMesh();
225
226 // Assign space for refractive indices cache and stripe effective indices
227 rsize = mesh->axis[0]->size();
228 zsize = mesh->axis[1]->size() + 1;
229 zbegin = 0;
230
231 if (geometry->isExtended(Geometry::DIRECTION_VERT, false) &&
232 abs(mesh->axis[1]->at(0) - geometry->getChild()->getBoundingBox().lower.c1) < SMALL)
233 zbegin = 1;
234 if (geometry->isExtended(Geometry::DIRECTION_TRAN, true) &&
235 abs(mesh->axis[0]->at(mesh->axis[0]->size() - 1) - geometry->getChild()->getBoundingBox().upper.c0) < SMALL)
236 --rsize;
237 if (geometry->isExtended(Geometry::DIRECTION_VERT, true) &&
238 abs(mesh->axis[1]->at(mesh->axis[1]->size() - 1) - geometry->getChild()->getBoundingBox().upper.c1) < SMALL)
239 --zsize;
240
241 nrCache.assign(rsize, std::vector<dcomplex, aligned_allocator<dcomplex>>(zsize));
242 ngCache.assign(rsize, std::vector<dcomplex, aligned_allocator<dcomplex>>(zsize));
243 veffs.resize(rsize);
244 nng.resize(rsize);
245
246 zfields.resize(zsize);
247
248 need_gain = false;
249 cache_outdated = true;
250 have_veffs = false;
251}
252
253 void EffectiveFrequencyCyl::onInvalidate() {
254 if (!modes.empty()) {
255 writelog(LOG_DETAIL, "Clearing computed modes");
256 modes.clear();
257 outWavelength.fireChanged();
258 outLoss.fireChanged();
259 outLightMagnitude.fireChanged();
260 outLightE.fireChanged();
261 }
262}
263
264/********* Here are the computations *********/
265
266 void EffectiveFrequencyCyl::updateCache() {
267 bool fresh = initCalculation();
268
269 if (fresh || cache_outdated || k0 != old_k0) {
270 // we need to update something
271
272 // Some additional checks
273 for (auto x : *mesh->axis[0]) {
274 if (x < 0.) throw BadMesh(getId(), "for cylindrical geometry no radial points can be negative");
275 }
276 if (abs(mesh->axis[0]->at(0)) > SMALL) throw BadMesh(getId(), "radial mesh must start from zero");
277
278 if (!modes.empty()) writelog(LOG_DETAIL, "Clearing computed modes");
279 modes.clear();
280
281 old_k0 = k0;
282
283 double lam = real(2e3 * PI / k0);
284
285 writelog(LOG_DEBUG, "Updating refractive indices cache");
286
287 double h = 1e6 * sqrt(SMALL);
288 double lam1 = lam - h, lam2 = lam + h;
289 double i2h = 0.5 / h;
290
291 shared_ptr<OrderedAxis> axis0, axis1;
292 {
293 shared_ptr<RectangularMesh<2>> midmesh = mesh->getElementMesh();
294 axis0 = plask::make_shared<OrderedAxis>(*midmesh->axis[0]);
295 axis1 = plask::make_shared<OrderedAxis>(*midmesh->axis[1]);
296 }
297 if (rsize == mesh->axis[0]->size())
298 axis0->addPoint(mesh->axis[0]->at(mesh->axis[0]->size() - 1) + 2. * OrderedAxis::MIN_DISTANCE);
299 if (zbegin == 0) axis1->addPoint(mesh->axis[1]->at(0) - 2. * OrderedAxis::MIN_DISTANCE);
300 if (zsize == mesh->axis[1]->size() + 1)
301 axis1->addPoint(mesh->axis[1]->at(mesh->axis[1]->size() - 1) + 2. * OrderedAxis::MIN_DISTANCE);
302
303 auto midmesh = plask::make_shared<RectangularMesh<2>>(axis0, axis1, mesh->getIterationOrder());
304 auto temp = inTemperature.hasProvider() ? inTemperature(midmesh) : LazyData<double>(midmesh->size(), 300.);
305 bool have_gain = false;
306 LazyData<Tensor2<double>> gain1, gain2;
307 auto carriers = inCarriersConcentration.hasProvider() ? inCarriersConcentration(CarriersConcentration::MAJORITY, midmesh)
308 : LazyData<double>(midmesh->size(), 0.);
309
310 for (size_t ir = 0; ir != rsize; ++ir) {
311 for (size_t iz = zbegin; iz < zsize; ++iz) {
312 size_t idx = midmesh->index(ir, iz - zbegin);
313 double T = temp[idx];
314 double cc = carriers[idx];
315 auto point = midmesh->at(idx);
316 auto material = geometry->getMaterial(point);
317 auto roles = geometry->getRolesAt(point);
318 // Nr = nr + i/(4π) λ g
319 // Ng = Nr - λ dN/dλ = Nr - λ dn/dλ - i/(4π) λ^2 dg/dλ
320 if (roles.find("QW") == roles.end() && roles.find("QD") == roles.end() && roles.find("gain") == roles.end()) {
321 nrCache[ir][iz] = material->Nr(lam, T, cc);
322 ngCache[ir][iz] = nrCache[ir][iz] - lam * (material->Nr(lam2, T, cc) - material->Nr(lam1, T, cc)) * i2h;
323 } else { // we ignore the material absorption as it should be considered in the gain already
324 need_gain = true;
325 if (!have_gain) {
326 gain1 = inGain(midmesh, lam1);
327 gain2 = inGain(midmesh, lam2);
328 have_gain = true;
329 }
330 double g = 0.5 * (gain1[idx].c00 + gain2[idx].c00);
331 double gs = (gain2[idx].c00 - gain1[idx].c00) * i2h;
332 double nr = real(material->Nr(lam, T, cc));
333 double ng = real(nr - lam * (material->Nr(lam2, T, cc) - material->Nr(lam1, T, cc)) * i2h);
334 nrCache[ir][iz] = dcomplex(nr, (0.25e-7 / PI) * lam * g);
335 ngCache[ir][iz] = dcomplex(ng, isnan(gs) ? 0. : -(0.25e-7 / PI) * lam * lam * gs);
336 }
337 }
338 if (zbegin != 0) {
339 nrCache[ir][0] = nrCache[ir][1];
340 ngCache[ir][0] = ngCache[ir][1];
341 }
342 }
343 cache_outdated = false;
344 have_veffs = false;
345 }
346}
347
348 void EffectiveFrequencyCyl::stageOne() {
349 updateCache();
350
351 if (!have_veffs) {
352 if (rstripe < 0) {
353 size_t main_stripe = getMainStripe();
354 // Compute effective frequencies for all stripes
355 std::exception_ptr error; // needed to handle exceptions from OMP loop
356 PLASK_OMP_PARALLEL_FOR
357 for (plask::openmp_size_t i = 0; i < rsize; ++i) {
358 if (error != std::exception_ptr()) continue; // just skip loops after error
359 try {
360 writelog(LOG_DETAIL, "Computing effective frequency for vertical stripe {0}", i);
361#ifndef NDEBUG
362 std::stringstream nrgs;
363 for (auto nr = nrCache[i].end(), ng = ngCache[i].end(); nr != nrCache[i].begin();) {
364 --nr;
365 --ng;
366 nrgs << ", (" << str(*nr) << ")/(" << str(*ng) << ")";
367 }
368 writelog(LOG_DEBUG, "Nr/Ng[{0}] = [{1} ]", i, nrgs.str().substr(1));
369#endif
370 dcomplex same_nr = nrCache[i].front();
371 dcomplex same_ng = ngCache[i].front();
372 bool all_the_same = true;
373 for (auto nr = nrCache[i].begin(), ng = ngCache[i].begin(); nr != nrCache[i].end(); ++nr, ++ng)
374 if (*nr != same_nr || *ng != same_ng) {
375 all_the_same = false;
376 break;
377 }
378 if (all_the_same) {
379 veffs[i] = 1.; // TODO make sure this is so!
380 nng[i] = same_nr * same_ng;
381 } else {
382 DataLog<dcomplex, dcomplex> log_stripe(getId(), format("stripe[{}]", i), "vlam", "det");
383 auto rootdigger = RootDigger::get(
384 this, [&](const dcomplex& x) { return this->detS1(2. - 4e3 * PI / x / k0, nrCache[i], ngCache[i]); },
385 log_stripe, stripe_root);
386 dcomplex start = (vlam == 0.) ? 2e3 * PI / k0 : vlam;
387 veffs[i] = freqv(rootdigger->find(start));
388 computeStripeNNg(i, i == main_stripe);
389 }
390 } catch (...) {
391#pragma omp critical
392 error = std::current_exception();
393 }
394 }
395 if (error != std::exception_ptr()) std::rethrow_exception(error);
396 } else {
397 // Compute effective frequencies for just one stripe
398 writelog(LOG_DETAIL, "Computing effective frequency for vertical stripe {0}", rstripe);
399#ifndef NDEBUG
400 std::stringstream nrgs;
401 for (auto nr = nrCache[rstripe].end(), ng = ngCache[rstripe].end(); nr != nrCache[rstripe].begin();) {
402 --nr;
403 --ng;
404 nrgs << ", (" << str(*nr) << ")/(" << str(*ng) << ")";
405 }
406 writelog(LOG_DEBUG, "Nr/Ng[{0}] = [{1} ]", rstripe, nrgs.str().substr(1));
407#endif
408 DataLog<dcomplex, dcomplex> log_stripe(getId(), format("stripe[{}]", rstripe), "vlam", "det");
409 auto rootdigger = RootDigger::get(
410 this, [&](const dcomplex& x) { return this->detS1(2. - 4e3 * PI / x / k0, nrCache[rstripe], ngCache[rstripe]); },
411 log_stripe, stripe_root);
412 dcomplex start = (vlam == 0.) ? 2e3 * PI / k0 : vlam;
413 veffs[rstripe] = freqv(rootdigger->find(start));
414 // Compute veffs and neffs for other stripes
415 computeStripeNNg(rstripe, true);
416 for (std::size_t i = 0; i < rsize; ++i)
417 if (i != std::size_t(rstripe)) computeStripeNNg(i);
418 }
419 assert(zintegrals.size() == zsize);
420
421#ifndef NDEBUG
422 std::stringstream strv;
423 for (size_t i = 0; i < veffs.size(); ++i) strv << ", " << str(veffs[i]);
424 writelog(LOG_DEBUG, "stripes veffs = [{0} ]", strv.str().substr(1));
425 std::stringstream strn;
426 for (size_t i = 0; i < nng.size(); ++i) strn << ", " << str(nng[i]);
427 writelog(LOG_DEBUG, "stripes <nng> = [{0} ]", strn.str().substr(1));
428#endif
429
430 have_veffs = true;
431
432 double rmin = INFINITY, rmax = -INFINITY, imin = INFINITY, imax = -INFINITY;
433 for (auto v : veffs) {
434 dcomplex lam = 2e3 * PI / (k0 * (1. - v / 2.));
435 if (real(lam) < rmin) rmin = real(lam);
436 if (real(lam) > rmax) rmax = real(lam);
437 if (imag(lam) < imin) imin = imag(lam);
438 if (imag(lam) > imax) imax = imag(lam);
439 }
440 writelog(LOG_DETAIL, "Wavelengths should be between {0}nm and {1}nm", str(dcomplex(rmin, imin)), str(dcomplex(rmax, imax)));
441 }
442}
443
444dcomplex EffectiveFrequencyCyl::detS1(const dcomplex& v,
445 const std::vector<dcomplex, aligned_allocator<dcomplex>>& NR,
446 const std::vector<dcomplex, aligned_allocator<dcomplex>>& NG,
447 std::vector<FieldZ>* saveto) {
448 if (saveto) (*saveto)[zbegin] = FieldZ(0., 1.);
449
450 std::vector<dcomplex> kz(zsize);
451 for (size_t i = zbegin; i < zsize; ++i) {
452 kz[i] = k0 * sqrt(NR[i] * NR[i] - v * NR[i] * NG[i]);
453 if (real(kz[i]) < 0.) kz[i] = -kz[i];
454 }
455
456 // dcomplex s1 = 1., s2 = 0., s3 = 0., s4 = 1.; // matrix S
457 //
458 // dcomplex phas = 1.;
459 // if (zbegin != 0)
460 // phas = exp(I * kz[zbegin] * (mesh->axis[1]->at(zbegin)-mesh->axis[1]->at(zbegin-1)));
461 //
462 // for (size_t i = zbegin+1; i < zsize; ++i) {
463 // // Compute shift inside one layer
464 // s1 *= phas;
465 // s3 *= phas * phas;
466 // s4 *= phas;
467 // // Compute matrix after boundary
468 // dcomplex p = 0.5 + 0.5 * kz[i] / kz[i-1];
469 // dcomplex m = 1.0 - p;
470 // dcomplex chi = 1. / (p - m * s3);
471 // // F0 = [ (-m*m + p*p)*chi*s1 m*s1*s4*chi + s2 ] [ F2 ]
472 // // B2 = [ (-m + p*s3)*chi s4*chi ] [ B0 ]
473 // s2 += s1*m*chi*s4;
474 // s1 *= (p*p - m*m) * chi;
475 // s3 = (p*s3-m) * chi;
476 // s4 *= chi;
477 // // Compute phase shift for the next step
478 // if (i != mesh->axis[1]->size())
479 // phas = exp(I * kz[i] * (mesh->axis[1]->at(i)-mesh->axis[1]->at(i-1)));
480 //
481 // // Compute fields
482 // if (saveto) {
483 // dcomplex F = -s2/s1, B = (s1*s4-s2*s3)/s1; // Assume F0 = 0 B0 = 1
484 // double aF = abs(F), aB = abs(B);
485 // // zero very small fields to avoid errors in plotting for long layers
486 // if (aF < 1e-8 * aB) F = 0.;
487 // if (aB < 1e-8 * aF) B = 0.;
488 // (*saveto)[i] = FieldZ(F, B);
489 // }
490 // }
491
492 MatrixZ T = MatrixZ::eye();
493 for (size_t i = zbegin; i < zsize - 1; ++i) {
494 double d;
495 if (i != zbegin || zbegin != 0)
496 d = mesh->axis[1]->at(i) - mesh->axis[1]->at(i - 1);
497 else
498 d = 0.;
499 dcomplex phas = exp(-I * kz[i] * d);
500 // Transfer through boundary
501 dcomplex n = 0.5 * kz[i] / kz[i + 1];
502 MatrixZ T1 = MatrixZ((0.5 + n), (0.5 - n), (0.5 - n), (0.5 + n));
503 T1.ff *= phas;
504 T1.fb /= phas;
505 T1.bf *= phas;
506 T1.bb /= phas;
507 T = T1 * T;
508 if (saveto) {
509 dcomplex F = T.fb, B = T.bb; // Assume F0 = 0 B0 = 1
510 double aF = abs(F), aB = abs(B);
511 // zero very small fields to avoid errors in plotting for long layers
512 if (aF < 1e-8 * aB) F = 0.;
513 if (aB < 1e-8 * aF) B = 0.;
514 (*saveto)[i + 1] = FieldZ(F, B);
515 }
516 }
517
518 if (saveto) {
519 dcomplex f;
520 if (emission == TOP) {
521 f = 1. / (*saveto)[zsize - 1].F / sqrt(NG[zsize - 1]);
522 (*saveto)[zsize - 1] = FieldZ(1., 0.);
523 } else {
524 // we dont have to scale fields as we have already set (*saveto)[zbegin] = FieldZ(0., 1.)
525 f = 1. / sqrt(NG[zbegin]);
526 (*saveto)[zsize - 1].B = 0.;
527 }
528 for (size_t i = zbegin; i < zsize - 1; ++i) (*saveto)[i] *= f;
529#ifndef NDEBUG
530 std::stringstream nrs;
531 for (size_t i = zbegin; i < zsize; ++i) nrs << "), (" << str((*saveto)[i].F) << ":" << str((*saveto)[i].B);
532 writelog(LOG_DEBUG, "vertical fields = [{0}) ]", nrs.str().substr(2));
533#endif
534 }
535
536 // return s4 - s2*s3/s1;
537
538 return T.bb; // F0 = 0 Bn = 0
539}
540
541 void EffectiveFrequencyCyl::computeStripeNNg(size_t stripe, bool save_integrals) {
542 size_t stripe0 = (rstripe < 0) ? stripe : rstripe;
543
544 nng[stripe] = 0.;
545
546 if (stripe != stripe0) veffs[stripe] = 0.;
547
548 std::vector<FieldZ> zfield(zsize);
549
550 dcomplex veff = veffs[stripe0];
551
552 // Compute fields
553 detS1(veff, nrCache[stripe0], ngCache[stripe0], &zfield);
554
555 double sum = 0.;
556
557 if (save_integrals) {
558#pragma omp critical
559 zintegrals.resize(zsize);
560 }
561
562 for (size_t i = zbegin + 1; i < zsize - 1; ++i) {
563 double d = mesh->axis[1]->at(i) - mesh->axis[1]->at(i - 1);
564 double weight = 0.;
565 dcomplex kz = k0 * sqrt(nrCache[stripe0][i] * nrCache[stripe0][i] - veff * nrCache[stripe0][i] * ngCache[stripe0][i]);
566 if (real(kz) < 0.) kz = -kz;
567 dcomplex w_ff, w_bb, w_fb, w_bf;
568 if (d != 0.) {
569 if (abs(imag(kz)) > SMALL) {
570 dcomplex kk = kz - conj(kz);
571 w_ff = (exp(-I * d * kk) - 1.) / kk;
572 w_bb = -(exp(+I * d * kk) - 1.) / kk;
573 } else
574 w_ff = w_bb = dcomplex(0., -d);
575 if (abs(real(kz)) > SMALL) {
576 dcomplex kk = kz + conj(kz);
577 w_fb = (exp(-I * d * kk) - 1.) / kk;
578 w_bf = -(exp(+I * d * kk) - 1.) / kk;
579 } else
580 w_ff = w_bb = dcomplex(0., -d);
581 weight = -imag(zfield[i].F * conj(zfield[i].F) * w_ff + zfield[i].F * conj(zfield[i].B) * w_fb +
582 zfield[i].B * conj(zfield[i].F) * w_bf + zfield[i].B * conj(zfield[i].B) * w_bb);
583 }
584 sum += weight;
585 if (save_integrals) {
586#pragma omp critical
587 zintegrals[i] = weight;
588 }
589 nng[stripe] += weight * nrCache[stripe][i] * ngCache[stripe][i];
590 if (stripe != stripe0) {
591 veffs[stripe] += weight * (nrCache[stripe][i] * nrCache[stripe][i] - nrCache[stripe0][i] * nrCache[stripe0][i]);
592 }
593 }
594
595 if (stripe != stripe0) {
596 veffs[stripe] += veff * nng[stripe0] * sum;
597 veffs[stripe] /= nng[stripe];
598 }
599
600 nng[stripe] /= sum;
601}
602
603 double EffectiveFrequencyCyl::integrateBessel(Mode& mode) {
604 double sum = 0;
605 for (size_t i = 0; i != rsize; ++i) {
606 double start = mesh->axis[0]->at(i);
607 double end = (i != rsize - 1) ? mesh->axis[0]->at(i + 1) : 3.0 * mesh->axis[0]->at(mesh->axis[0]->size() - 1);
608 double err = perr;
609 mode.rweights[i] = patterson<double, double>([this, &mode](double r) { return r * abs2(mode.rField(r)); }, start, end, err);
610 // TODO use exponential asymptotic approximation to compute weight in the last stripe
611 sum += mode.rweights[i];
612 }
613 // TODO consider m <> 0
614 double f = 1e12 / sum;
615 for (double& w : mode.rweights) w *= f;
616 return 2. * PI * sum;
617}
618
619 void EffectiveFrequencyCyl::computeBessel(size_t i,
620 dcomplex v,
621 const Mode& mode,
622 dcomplex* J1,
623 dcomplex* H1,
624 dcomplex* J2,
625 dcomplex* H2) {
626 double r = mesh->axis[0]->at(i);
627 dcomplex x1 = r * k0 * sqrt(nng[i - 1] * (veffs[i - 1] - v));
628 if (real(x1) < 0.) x1 = -x1;
629 if (imag(x1) > SMALL) x1 = -x1;
630 dcomplex x2 = r * k0 * sqrt(nng[i] * (veffs[i] - v));
631 if (real(x2) < 0.) x2 = -x2;
632 if (imag(x2) > SMALL) x2 = -x2;
633 // Compute Bessel functions and their derivatives
634 double Jr[2], Ji[2], Hr[2], Hi[2];
635 long nz, ierr;
636 zbesj(x1.real(), x1.imag(), mode.m, 1, 2, Jr, Ji, nz, ierr);
637 if (ierr != 0)
638 throw ComputationError(getId(), "could not compute J({0}, {1})\n @ r = {2} um, lam = {3} nm, vlam = {4} nm", mode.m,
639 str(x1), r, str(lambda(v)), str(lambda(veffs[i - 1])));
640 zbesh(x1.real(), x1.imag(), mode.m, 1, MH, 2, Hr, Hi, nz, ierr);
641 if (ierr != 0)
642 throw ComputationError(getId(), "could not compute H({0}, {1})\n @ r = {2} um, lam = {3} nm, vlam = {4} nm", mode.m,
643 str(x1), r, str(lambda(v)), str(lambda(veffs[i - 1])));
644 for (int j = 0; j < 2; ++j) {
645 J1[j] = dcomplex(Jr[j], Ji[j]);
646 H1[j] = dcomplex(Hr[j], Hi[j]);
647 }
648 zbesj(x2.real(), x2.imag(), mode.m, 1, 2, Jr, Ji, nz, ierr);
649 if (ierr != 0)
650 throw ComputationError(getId(), "could not compute J({0}, {1})\n @ r = {2} um, lam = {3} nm, vlam = {4} nm", mode.m,
651 str(x1), r, str(lambda(v)), str(lambda(veffs[i])));
652 zbesh(x2.real(), x2.imag(), mode.m, 1, MH, 2, Hr, Hi, nz, ierr);
653 if (ierr != 0)
654 throw ComputationError(getId(), "could not compute H({0}, {1})\n @ r = {2} um, lam = {3} nm, vlam = {4} nm", mode.m,
655 str(x1), r, str(lambda(v)), str(lambda(veffs[i])));
656 for (int j = 0; j < 2; ++j) {
657 J2[j] = dcomplex(Jr[j], Ji[j]);
658 H2[j] = dcomplex(Hr[j], Hi[j]);
659 }
660 J1[1] *= x1;
661 H1[1] *= x1;
662 J2[1] *= x2;
663 H2[1] *= x2;
664}
665
666#define zgeev F77_GLOBAL(zgeev, ZGEEV)
667 F77SUB zgeev(const char& jobvl,
668 const char& jobvr,
669 const int& n,
670 dcomplex* a,
671 const int& lda,
672 dcomplex* w,
673 dcomplex* vl,
674 const int& ldvl,
675 dcomplex* vr,
676 const int& ldvr,
677 dcomplex* work,
678 const int& lwork,
679 double* rwork,
680 int& info);
681
682dcomplex EffectiveFrequencyCyl::detS(const dcomplex& lam, Mode& mode, bool save) {
683 dcomplex v = freqv(lam);
684 dcomplex J1[2], H1[2];
685 dcomplex J2[2], H2[2];
686
687 double m = mode.m;
688
689 if (determinant == DETERMINANT_FULL) {
690 const size_t N = 2 * rsize;
691 if (!save) {
692 ZgbMatrix matrix(N);
693 // In the innermost layer the solution is only the J function
694 matrix(0, 0) = 1.;
695 matrix(0, 1) = 0.;
696 for (size_t i = 1, n = 1; i != rsize; ++i, n += 2) {
697 computeBessel(i, v, mode, J1, H1, J2, H2);
698 matrix(n - 1, n + 1) = 0.;
699 matrix(n, n - 1) = H1[0];
700 matrix(n, n) = J1[0];
701 matrix(n, n + 1) = -H2[0];
702 matrix(n, n + 2) = -J2[0];
703 matrix(n + 1, n - 1) = m * H1[0] - H1[1];
704 matrix(n + 1, n) = m * J1[0] - J1[1];
705 matrix(n + 1, n + 1) = -m * H2[0] + H2[1];
706 matrix(n + 1, n + 2) = -m * J2[0] + J2[1];
707 matrix(n + 2, n) = 0.;
708 }
709 // In the outermost area there must only outgoing wave, so J = 0.
710 matrix(N - 1, N - 2) = 0.;
711 matrix(N - 1, N - 1) = 1.;
712 return matrix.determinant();
713 } else {
714 const std::size_t NN = N * N;
715 dcomplex *matrix = aligned_malloc<dcomplex>(NN), *evals = aligned_malloc<dcomplex>(N),
716 *evecs = aligned_malloc<dcomplex>(NN);
717 aligned_unique_ptr<dcomplex[]> pmatrix(matrix), pevals(evals), pevecs(evecs);
718
719 std::fill_n(matrix, N * N, 0.);
720 // In the innermost layer the solution is only the J function
721 matrix[0] = 1.;
722 matrix[N] = 0.;
723 for (size_t i = 1, n = 1; i != rsize; ++i, n += 2) {
724 computeBessel(i, v, mode, J1, H1, J2, H2);
725 matrix[n + N * (n - 1)] = H1[0];
726 matrix[n + N * (n)] = J1[0];
727 matrix[n + N * (n + 1)] = -H2[0];
728 matrix[n + N * (n + 2)] = -J2[0];
729 matrix[n + 1 + N * (n - 1)] = m * H1[0] - H1[1];
730 matrix[n + 1 + N * (n)] = m * J1[0] - J1[1];
731 matrix[n + 1 + N * (n + 1)] = -m * H2[0] + H2[1];
732 matrix[n + 1 + N * (n + 2)] = -m * J2[0] + J2[1];
733 }
734 // In the outermost area there must only outgoing wave, so J = 0.
735 matrix[NN - 1] = 1.;
736 const std::size_t lwork = 2 * N + 1;
737 aligned_unique_ptr<dcomplex[]> work(aligned_malloc<dcomplex>(lwork));
738 aligned_unique_ptr<double[]> rwork(aligned_malloc<double>(2 * N));
739 int info;
740 zgeev('N', 'V', int(N), matrix, int(N), evals, nullptr, 1, evecs, int(N), work.get(), int(lwork), rwork.get(), info);
741 if (info != 0) throw ComputationError(getId(), "could not compute eigenvalues of radial continuity matrix");
742 // Find the eigenvalue with the smallest absolute value
743 double minabs = INFINITY;
744 size_t minidx;
745 dcomplex res = 1.;
746 for (int i = 0; i != N; ++i) {
747 res *= evals[i];
748 double abs = abs2(evals[i]);
749 if (abs < minabs) {
750 minabs = abs;
751 minidx = i;
752 }
753 }
754 dcomplex* evec = evecs + N * minidx;
755 for (size_t i = 0, n = 0; i < rsize; ++i, n += 2) {
756 mode.rfields[i] = FieldR(evec[n + 1], evec[n]);
757 }
758 dcomplex f = 1e6 * sqrt(1. / integrateBessel(mode)); // 1e6: V/µm -> V/m
759 for (size_t r = 0; r != rsize; ++r) mode.rfields[r] *= f;
760 return res;
761 }
762 } else if (determinant == DETERMINANT_INWARDS) {
763 MatrixR T = MatrixR::eye();
764 mode.rfields[rsize - 1] = FieldR(0., 1.);
765 for (size_t i = rsize - 1; i > 0; --i) {
766 computeBessel(i, v, mode, J1, H1, J2, H2);
767 MatrixR M1(J1[0], H1[0], m * J1[0] - J1[1], m * H1[0] - H1[1]);
768 MatrixR M2(J2[0], H2[0], m * J2[0] - J2[1], m * H2[0] - H2[1]);
769 T = M1.solve(M2) * T;
770 if (save) mode.rfields[i - 1] = T * mode.rfields[rsize - 1];
771 }
772 if (save) {
773 dcomplex f = 1e6 * sqrt(1. / integrateBessel(mode)); // 1e6: V/µm -> V/m
774 for (size_t r = 0; r != rsize; ++r) mode.rfields[r] *= f;
775 }
776 // In the innermost area there must not be any infinity, so H_0 = 0.
777 // j = [ JJ JH ] 0
778 // 0 = [ HJ HH ] 1
779 return T.HH / T.JH * H1[0] / J1[0];
780 } else {
781 MatrixR T = MatrixR::eye();
782 mode.rfields[0] = FieldR(1., 0.);
783 for (size_t i = 1; i < rsize; ++i) {
784 computeBessel(i, v, mode, J1, H1, J2, H2);
785 MatrixR M1(J1[0], H1[0], m * J1[0] - J1[1], m * H1[0] - H1[1]);
786 MatrixR M2(J2[0], H2[0], m * J2[0] - J2[1], m * H2[0] - H2[1]);
787 T = M2.solve(M1) * T;
788 if (save) mode.rfields[i] = T * mode.rfields[0];
789 }
790 if (save) {
791 dcomplex f = 1e6 * sqrt(1. / integrateBessel(mode)); // 1e6: V/µm -> V/m
792 for (size_t r = 0; r != rsize; ++r) mode.rfields[r] *= f;
793 }
794 // In the outermost area there is only outgoing wave, so J_N = 0.
795 // 0 = [ JJ JH ] 1
796 // h = [ HJ HH ] 0
797 return T.JJ / T.HJ * J2[0] / H2[0];
798 }
799}
800
801 double EffectiveFrequencyCyl::getTotalAbsorption(Mode& mode) {
802 if (!mode.have_fields) {
803 size_t stripe = getMainStripe();
804 detS1(veffs[stripe], nrCache[stripe], ngCache[stripe], &zfields); // compute vertical part
805 detS(mode.lam, mode, true); // compute horizontal part
806 mode.have_fields = true;
807 }
808
809 double result = 0.;
810 dcomplex lam0 = 2e3 * PI / k0;
811
812 for (size_t ir = 0; ir < rsize; ++ir) {
813 for (size_t iz = zbegin + 1; iz < zsize - 1; ++iz) {
814 dcomplex n = nrCache[ir][iz] + ngCache[ir][iz] * (1. - mode.lam / lam0);
815 double absp = -2. * real(n) * imag(n);
816 result += absp * mode.rweights[ir] * zintegrals[iz]; // [dV] = µm³
817 // double err = 1e-6, erz = 1e-6;
818 // double rstart = mesh->axis[0][ir];
819 // double rend = (ir != rsize-1)? mesh->axis[0][ir+1] : 3.0 * mesh->axis[0][ir];
820 // result += absp * 2.*PI *
821 // patterson<double>([this,&mode](double r){return r * abs2(mode.rField(r));}, rstart, rend, err) *
822 // patterson<double>([&](double z){
823 // size_t stripe = getMainStripe();
824 // dcomplex kz = k0 * sqrt(nrCache[stripe][iz]*nrCache[stripe][iz] - veffs[stripe] *
825 // nrCache[stripe][iz]*ngCache[stripe][iz]); if (real(kz) < 0.) kz = -kz; z -= mesh->axis[1][iz]; dcomplex phasz
826 // = exp(- I * kz * z); return abs2(zfields[iz].F * phasz + zfields[iz].B / phasz);
827 // }, mesh->axis[1][iz-1], mesh->axis[1][iz], erz);
828 }
829 }
830 result *= 2e-9 * PI / real(mode.lam) * mode.power; // 1e-9: µm³ / nm -> m², 2: ½ is already hidden in mode.power
831 return result;
832}
833
834 double EffectiveFrequencyCyl::getTotalAbsorption(size_t num) {
835 if (modes.size() <= num || k0 != old_k0) throw NoValue("absorption");
836
837 return getTotalAbsorption(modes[num]);
838}
839
840 double EffectiveFrequencyCyl::getGainIntegral(Mode& mode) {
841 if (!mode.have_fields) {
842 size_t stripe = getMainStripe();
843 detS1(veffs[stripe], nrCache[stripe], ngCache[stripe], &zfields); // compute vertical part
844 detS(mode.lam, mode, true); // compute horizontal part
845 mode.have_fields = true;
846 }
847
848 double result = 0.;
849 dcomplex lam0 = 2e3 * PI / k0;
850
851 auto midmesh = mesh->getElementMesh();
852
853 for (size_t ir = 0; ir < rsize; ++ir) {
854 for (size_t iz = zbegin + 1; iz < zsize - 1; ++iz) {
855 auto roles = geometry->getRolesAt(midmesh->at(ir, iz - 1));
856 if (roles.find("QW") != roles.end() || roles.find("QD") != roles.end() || roles.find("gain") != roles.end()) {
857 dcomplex n = nrCache[ir][iz] + ngCache[ir][iz] * (1. - mode.lam / lam0);
858 double absp = -2. * real(n) * imag(n);
859 result += absp * mode.rweights[ir] * zintegrals[iz]; // [dV] = µm³
860 }
861 }
862 }
863 result *= 2e-9 * PI / real(mode.lam) * mode.power; // 1e-9: µm³ / nm -> m², 2: ½ is already hidden in mode.power
864 return -result;
865}
866
867 double EffectiveFrequencyCyl::getGainIntegral(size_t num) {
868 if (modes.size() <= num || k0 != old_k0) throw NoValue("absorption");
869
870 return getGainIntegral(modes[num]);
871}
872
873 template <typename FieldT> struct EffectiveFrequencyCyl::FieldDataBase : public LazyDataImpl<FieldT> {
874 EffectiveFrequencyCyl* solver;
875 std::size_t num;
876
877 FieldDataBase(EffectiveFrequencyCyl* solver, std::size_t num);
878
879 protected:
880 inline FieldT value(dcomplex val) const;
881 double scale;
882};
883
884 template <>
885 EffectiveFrequencyCyl::FieldDataBase<double>::FieldDataBase(EffectiveFrequencyCyl* solver, std::size_t num)
886 : solver(solver), num(num), scale(1e-3 * solver->modes[num].power) {}
887
888 template <> double EffectiveFrequencyCyl::FieldDataBase<double>::value(dcomplex val) const { return scale * abs2(val); }
889
890 template <>
891 EffectiveFrequencyCyl::FieldDataBase<Vec<3, dcomplex>>::FieldDataBase(EffectiveFrequencyCyl* solver, std::size_t num)
892 : solver(solver),
893 num(num),
894 scale(sqrt(2e-3 * phys::Z0 * solver->modes[num].power))
895// <M> = ½ E conj(E) / Z0
896{}
897
898 template <> Vec<3, dcomplex> EffectiveFrequencyCyl::FieldDataBase<Vec<3, dcomplex>>::value(dcomplex val) const {
899 return Vec<3, dcomplex>(0., scale * val, 0.);
900}
901
902 template <typename FieldT> struct EffectiveFrequencyCyl::FieldDataInefficient : public FieldDataBase<FieldT> {
903 shared_ptr<const MeshD<2>> dst_mesh;
904 size_t stripe;
905
906 FieldDataInefficient(EffectiveFrequencyCyl* solver, std::size_t num, const shared_ptr<const MeshD<2>>& dst_mesh, size_t stripe)
907 : FieldDataBase<FieldT>(solver, num), dst_mesh(dst_mesh), stripe(stripe) {}
908
909 size_t size() const override { return dst_mesh->size(); }
910
911 FieldT at(size_t id) const override {
912 auto point = dst_mesh->at(id);
913 double r = point.c0;
914 double z = point.c1;
915 if (r < 0) r = -r;
916
917 dcomplex val = this->solver->modes[this->num].rField(r);
918
919 size_t iz = this->solver->mesh->axis[1]->findIndex(z);
920 if (iz >= this->solver->zsize)
921 iz = this->solver->zsize - 1;
922 else if (iz < this->solver->zbegin)
923 iz = this->solver->zbegin;
924 dcomplex kz = this->solver->k0 *
925 sqrt(this->solver->nrCache[stripe][iz] * this->solver->nrCache[stripe][iz] -
926 this->solver->veffs[stripe] * this->solver->nrCache[stripe][iz] * this->solver->ngCache[stripe][iz]);
927 if (real(kz) < 0.) kz = -kz;
928 z -= this->solver->mesh->axis[1]->at(max(int(iz) - 1, 0));
929 dcomplex phasz = exp(-I * kz * z);
930 val *= this->solver->zfields[iz].F * phasz + this->solver->zfields[iz].B / phasz;
931
932 return this->value(val);
933 }
934};
935
936 template <typename FieldT> struct EffectiveFrequencyCyl::FieldDataEfficient : public FieldDataBase<FieldT> {
937 shared_ptr<const RectangularMesh<2>> rect_mesh;
938 std::vector<dcomplex> valr, valz;
939
940 FieldDataEfficient(EffectiveFrequencyCyl* solver,
941 std::size_t num,
942 const shared_ptr<const RectangularMesh<2>>& rect_mesh,
943 size_t stripe)
944 : FieldDataBase<FieldT>(solver, num),
945 rect_mesh(rect_mesh),
946 valr(rect_mesh->axis[0]->size()),
947 valz(rect_mesh->axis[1]->size()) {
948 std::exception_ptr error; // needed to handle exceptions from OMP loop
949
950 PLASK_OMP_PARALLEL
951 {
952#pragma omp for nowait
953 for (int idr = 0; idr < int(rect_mesh->axis[0]->size());
954 ++idr) { // idr can't be size_t since MSVC does not support omp newer than 2
955 if (error) continue;
956 double r = rect_mesh->axis[0]->at(idr);
957 if (r < 0.) r = -r;
958 try {
959 valr[idr] = solver->modes[num].rField(r);
960 } catch (...) {
961#pragma omp critical
962 error = std::current_exception();
963 }
964 }
965
966 if (!error) {
967#pragma omp for
968 for (int idz = 0; idz < int(rect_mesh->axis[1]->size());
969 ++idz) { // idz can't be size_t since MSVC does not support omp newer than 2
970 double z = rect_mesh->axis[1]->at(idz);
971 size_t iz = solver->mesh->axis[1]->findIndex(z);
972 if (iz >= solver->zsize)
973 iz = solver->zsize - 1;
974 else if (iz < solver->zbegin)
975 iz = solver->zbegin;
976 dcomplex kz =
977 solver->k0 * sqrt(solver->nrCache[stripe][iz] * solver->nrCache[stripe][iz] -
978 solver->veffs[stripe] * solver->nrCache[stripe][iz] * solver->ngCache[stripe][iz]);
979 if (real(kz) < 0.) kz = -kz;
980 z -= solver->mesh->axis[1]->at(max(int(iz) - 1, 0));
981 dcomplex phasz = exp(-I * kz * z);
982 valz[idz] = solver->zfields[iz].F * phasz + solver->zfields[iz].B / phasz;
983 }
984 }
985 }
986 if (error) std::rethrow_exception(error);
987 }
988
989 size_t size() const override { return rect_mesh->size(); }
990
991 FieldT at(size_t idx) const override {
992 size_t i0 = rect_mesh->index0(idx);
993 size_t i1 = rect_mesh->index1(idx);
994 return this->value(valr[i0] * valz[i1]);
995 }
996
997 DataVector<const FieldT> getAll() const override {
998 DataVector<FieldT> results(rect_mesh->size());
999
1000 if (rect_mesh->getIterationOrder() == RectangularMesh<2>::ORDER_10) {
1001 PLASK_OMP_PARALLEL_FOR
1002 for (plask::openmp_size_t i1 = 0; i1 < rect_mesh->axis[1]->size(); ++i1) {
1003 FieldT* data = results.data() + i1 * rect_mesh->axis[0]->size();
1004 for (size_t i0 = 0; i0 < rect_mesh->axis[0]->size(); ++i0) {
1005 dcomplex f = valr[i0] * valz[i1];
1006 data[i0] = this->value(f);
1007 }
1008 }
1009 } else {
1010 PLASK_OMP_PARALLEL_FOR
1011 for (plask::openmp_size_t i0 = 0; i0 < rect_mesh->axis[0]->size(); ++i0) {
1012 FieldT* data = results.data() + i0 * rect_mesh->axis[1]->size();
1013 for (size_t i1 = 0; i1 < rect_mesh->axis[1]->size(); ++i1) {
1014 dcomplex f = valr[i0] * valz[i1];
1015 data[i1] = this->value(f);
1016 }
1017 }
1018 }
1019 return results;
1020 }
1021};
1022
1023 const LazyData<double> EffectiveFrequencyCyl::getLightMagnitude(std::size_t num,
1024 const shared_ptr<const MeshD<2>>& dst_mesh,
1025 InterpolationMethod) {
1026 this->writelog(LOG_DEBUG, "Getting light magnitude");
1027
1028 if (modes.size() <= num || k0 != old_k0) throw NoValue(LightMagnitude::NAME);
1029
1030 size_t stripe = getMainStripe();
1031
1032 if (!modes[num].have_fields) {
1033 // Compute vertical part
1034 detS1(veffs[stripe], nrCache[stripe], ngCache[stripe], &zfields);
1035 // Compute horizontal part
1036 detS(modes[num].lam, modes[num], true);
1037#ifndef NDEBUG
1038 {
1039 std::stringstream nrs;
1040 for (size_t i = 0; i < rsize; ++i) nrs << "), (" << str(modes[num].rfields[i].J) << ":" << str(modes[num].rfields[i].H);
1041 writelog(LOG_DEBUG, "horizontal fields = [{0}) ]", nrs.str().substr(2));
1042 }
1043#endif
1044 modes[num].have_fields = true;
1045 }
1046
1047 if (auto rect_mesh = dynamic_pointer_cast<const RectangularMesh<2>>(dst_mesh))
1048 return LazyData<double>(new FieldDataEfficient<double>(this, num, rect_mesh, stripe));
1049 else
1050 return LazyData<double>(new FieldDataInefficient<double>(this, num, dst_mesh, stripe));
1051}
1052
1053 const LazyData<Vec<3, dcomplex>> EffectiveFrequencyCyl::getElectricField(std::size_t num,
1054 const shared_ptr<const MeshD<2>>& dst_mesh,
1055 InterpolationMethod) {
1056 this->writelog(LOG_DEBUG, "Getting light electric field");
1057
1058 if (modes.size() <= num || k0 != old_k0) throw NoValue(LightMagnitude::NAME);
1059
1060 size_t stripe = getMainStripe();
1061
1062 if (!modes[num].have_fields) {
1063 // Compute vertical part
1064 detS1(veffs[stripe], nrCache[stripe], ngCache[stripe], &zfields);
1065 // Compute horizontal part
1066 detS(modes[num].lam, modes[num], true);
1067#ifndef NDEBUG
1068 {
1069 std::stringstream nrs;
1070 for (size_t i = 0; i < rsize; ++i) nrs << "), (" << str(modes[num].rfields[i].J) << ":" << str(modes[num].rfields[i].H);
1071 writelog(LOG_DEBUG, "horizontal fields = [{0}) ]", nrs.str().substr(2));
1072 }
1073#endif
1074 modes[num].have_fields = true;
1075 }
1076
1077 if (auto rect_mesh = dynamic_pointer_cast<const RectangularMesh<2>>(dst_mesh))
1078 return LazyData<Vec<3, dcomplex>>(new FieldDataEfficient<Vec<3, dcomplex>>(this, num, rect_mesh, stripe));
1079 else
1080 return LazyData<Vec<3, dcomplex>>(new FieldDataInefficient<Vec<3, dcomplex>>(this, num, dst_mesh, stripe));
1081}
1082
1083 const LazyData<dcomplex> EffectiveFrequencyCyl::getRefractiveIndex(RefractiveIndex::EnumType,
1084 const shared_ptr<const MeshD<2>>& dst_mesh,
1085 dcomplex lam,
1086 InterpolationMethod) {
1087 if (!isnan(real(lam))) throw BadInput(this->getId(), "wavelength cannot be specified for outRefractiveIndex in this solver");
1088 this->writelog(LOG_DEBUG, "Getting refractive indices");
1089 dcomplex lam0 = 2e3 * PI / k0;
1090 updateCache();
1091 InterpolationFlags flags(geometry);
1092 return LazyData<dcomplex>(dst_mesh->size(), [this, dst_mesh, flags, lam0](size_t j) -> dcomplex {
1093 auto point = flags.wrap(dst_mesh->at(j));
1094 size_t ir = this->mesh->axis[0]->findIndex(point.c0);
1095 if (ir != 0) --ir;
1096 if (ir >= this->rsize) ir = this->rsize - 1;
1097 size_t iz = this->mesh->axis[1]->findIndex(point.c1);
1098 if (iz < this->zbegin)
1099 iz = this->zbegin;
1100 else if (iz >= zsize)
1101 iz = this->zsize - 1;
1102 return this->nrCache[ir][iz] /* + this->ngCache[ir][iz] * (1. - lam/lam0)*/;
1103 });
1104}
1105
1106 struct EffectiveFrequencyCyl::HeatDataImpl : public LazyDataImpl<double> {
1107 EffectiveFrequencyCyl* solver;
1108 shared_ptr<const MeshD<2>> dest_mesh;
1109 InterpolationFlags flags;
1110 std::vector<LazyData<double>> EE;
1111 dcomplex lam0;
1112
1113 HeatDataImpl(EffectiveFrequencyCyl* solver, const shared_ptr<const MeshD<2>>& dst_mesh, InterpolationMethod method)
1114 : solver(solver), dest_mesh(dst_mesh), flags(solver->geometry), EE(solver->modes.size()), lam0(2e3 * PI / solver->k0) {
1115 for (size_t m = 0; m != solver->modes.size(); ++m) EE[m] = solver->getLightMagnitude(m, dst_mesh, method);
1116 }
1117
1118 size_t size() const override { return dest_mesh->size(); }
1119
1120 double at(size_t j) const override {
1121 double result = 0.;
1122 auto point = flags.wrap(dest_mesh->at(j));
1123 size_t ir = solver->mesh->axis[0]->findIndex(point.c0);
1124 if (ir != 0) --ir;
1125 if (ir >= solver->rsize) ir = solver->rsize - 1;
1126 size_t iz = solver->mesh->axis[1]->findIndex(point.c1);
1127 if (iz < solver->zbegin)
1128 iz = solver->zbegin;
1129 else if (iz >= solver->zsize)
1130 iz = solver->zsize - 1;
1131 for (size_t m = 0; m != solver->modes.size(); ++m) { // we sum heats from all modes
1132 dcomplex n = solver->nrCache[ir][iz] + solver->ngCache[ir][iz] * (1. - solver->modes[m].lam / lam0);
1133 double absp = -2. * real(n) * imag(n);
1134 result += 2e9 * PI / real(solver->modes[m].lam) * absp * EE[m][j]; // 1e9: 1/nm -> 1/m
1135 }
1136 return result;
1137 }
1138};
1139
1140 const LazyData<double> EffectiveFrequencyCyl::getHeat(const shared_ptr<const MeshD<2>>& dst_mesh, InterpolationMethod method) {
1141 // This is somehow naive implementation using the field value from the mesh points. The heat may be slightly off
1142 // in case of fast varying light intensity and too sparse mesh.
1143
1144 writelog(LOG_DEBUG, "Getting heat absorbed from {0} mode{1}", modes.size(), (modes.size() == 1) ? "" : "s");
1145 if (modes.size() == 0) return LazyData<double>(dst_mesh->size(), 0.);
1146 return LazyData<double>(new HeatDataImpl(this, dst_mesh, method));
1147}
1148
1149}}} // namespace plask::optical::effective