PLaSK library
Loading...
Searching...
No Matches
eim.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 "eim.hpp"
15#include <exception>
16
17namespace plask { namespace optical { namespace effective {
18
19EffectiveIndex2D::EffectiveIndex2D(const std::string& name)
20 : SolverWithMesh<Geometry2DCartesian, RectangularMesh<2>>(name),
21 log_value(dataLog<dcomplex, dcomplex>("Neff", "Neff", "det")),
22 stripex(0.),
23 polarization(TE),
24 recompute_neffs(true),
25 emission(FRONT),
26 vneff(0.),
27 outNeff(this, &EffectiveIndex2D::getEffectiveIndex, &EffectiveIndex2D::nmodes),
28 outLightMagnitude(this, &EffectiveIndex2D::getLightMagnitude, &EffectiveIndex2D::nmodes),
29 outLightE(this, &EffectiveIndex2D::getElectricField, &EffectiveIndex2D::nmodes),
30 outRefractiveIndex(this, &EffectiveIndex2D::getRefractiveIndex),
31 outHeat(this, &EffectiveIndex2D::getHeat),
32 k0(2e3 * PI / 980) {
33 inTemperature = 300.;
34 inGain = NAN;
35 root.tolx = 1.0e-6;
36 root.tolf_min = 1.0e-7;
37 root.tolf_max = 2.0e-5;
38 root.maxiter = 500;
39 root.method = RootDigger::ROOT_MULLER;
40 stripe_root.tolx = 1.0e-6;
41 stripe_root.tolf_min = 1.0e-7;
42 stripe_root.tolf_max = 1.0e-5;
43 stripe_root.maxiter = 500;
44 stripe_root.method = RootDigger::ROOT_MULLER;
45 inTemperature.changedConnectMethod(this, &EffectiveIndex2D::onInputChange);
46 inGain.changedConnectMethod(this, &EffectiveIndex2D::onInputChange);
47 inCarriersConcentration.changedConnectMethod(this, &EffectiveIndex2D::onInputChange);
48}
49
50void EffectiveIndex2D::loadConfiguration(XMLReader& reader, Manager& manager) {
51 while (reader.requireTagOrEnd()) {
52 std::string param = reader.getNodeName();
53 if (param == "mode") {
54 auto pol = reader.getAttribute("polarization");
55 if (pol) {
56 if (*pol == "TE")
57 polarization = TE;
58 else if (*pol == "TM")
59 polarization = TM;
60 else
61 throw BadInput(getId(), "wrong polarization specification '{0}' in XML", *pol);
62 }
63 k0 = 2e3 * PI / reader.getAttribute<double>("wavelength", real(2e3 * PI / k0));
64 stripex = reader.getAttribute<double>("vat", stripex);
65 vneff = reader.getAttribute<dcomplex>("vneff", vneff);
66 emission = reader.enumAttribute<Emission>("emission").value("front", FRONT).value("back", BACK).get(emission);
67 reader.requireTagEnd();
68 } else if (param == "root") {
69 RootDigger::readRootDiggerConfig(reader, root);
70 } else if (param == "stripe-root") {
71 RootDigger::readRootDiggerConfig(reader, stripe_root);
72 } else if (param == "mirrors") {
73 double R1 = reader.requireAttribute<double>("R1");
74 double R2 = reader.requireAttribute<double>("R2");
75 mirrors.reset(std::make_pair(R1, R2));
76 reader.requireTagEnd();
77 } else if (param == "mesh") {
78 auto name = reader.getAttribute("ref");
79 if (!name)
80 name.reset(reader.requireTextInCurrentTag());
81 else
82 reader.requireTagEnd();
83 auto found = manager.meshes.find(*name);
84 if (found != manager.meshes.end()) {
85 auto mesh1 = dynamic_pointer_cast<MeshAxis>(found->second);
86 auto mesh2 = dynamic_pointer_cast<RectangularMesh<2>>(found->second);
87 if (mesh1)
88 this->setHorizontalMesh(mesh1);
89 else if (mesh2)
90 this->setMesh(mesh2);
91 else {
92 auto generator1 = dynamic_pointer_cast<MeshGeneratorD<1>>(found->second);
93 auto generator2 = dynamic_pointer_cast<MeshGeneratorD<2>>(found->second);
94 if (generator1)
95 this->setMesh(plask::make_shared<RectangularMesh2DFrom1DGenerator>(generator1));
96 else if (generator2)
97 this->setMesh(generator2);
98 else
99 throw BadInput(this->getId(), "mesh or generator '{0}' of wrong type", *name);
100 }
101 }
102 } else
103 parseStandardConfiguration(reader, manager, "<geometry>, <mesh>, <mode>, <root>, <stripe-root>, or <outer>");
104 }
105}
106
107std::vector<dcomplex> EffectiveIndex2D::searchVNeffs(dcomplex neff1, dcomplex neff2, size_t resteps, size_t imsteps, dcomplex eps) {
108 updateCache();
109
110 size_t stripe = mesh->tran()->findIndex(stripex);
111 if (stripe < xbegin)
112 stripe = xbegin;
113 else if (stripe >= xend)
114 stripe = xend - 1;
115
116 if (eps.imag() == 0.) eps.imag(eps.real());
117
118 if (real(eps) <= 0. || imag(eps) <= 0.) throw BadInput(this->getId(), "bad precision specified");
119
120 double re0 = real(neff1), im0 = imag(neff1);
121 double re1 = real(neff2), im1 = imag(neff2);
122 if (re0 > re1) std::swap(re0, re1);
123 if (im0 > im1) std::swap(im0, im1);
124
125 if (re0 == 0. && re1 == 0.) {
126 re0 = 1e30;
127 re1 = -1e30;
128 for (size_t i = ybegin; i != yend; ++i) {
129 dcomplex n = nrCache[stripe][i];
130 if (n.real() < re0) re0 = n.real();
131 if (n.real() > re1) re1 = n.real();
132 }
133 } else if (re0 == 0. || re1 == 0.)
134 throw BadInput(getId(), "bad area to browse specified");
135 if (im0 == 0. && im1 == 0.) {
136 im0 = 1e30;
137 im1 = -1e30;
138 for (size_t i = ybegin; i != yend; ++i) {
139 dcomplex n = nrCache[stripe][i];
140 if (n.imag() < im0) im0 = n.imag();
141 if (n.imag() > im1) im1 = n.imag();
142 }
143 }
144 neff1 = dcomplex(re0, im0);
145 neff2 = dcomplex(re1, im1);
146
147 auto ranges = findZeros(
148 this, [&](const dcomplex& z) { return this->detS1(z, nrCache[stripe]); }, neff1, neff2, resteps, imsteps, eps);
149 std::vector<dcomplex> results;
150 results.reserve(ranges.size());
151 for (auto zz : ranges) results.push_back(0.5 * (zz.first + zz.second));
152
153 if (maxLoglevel >= LOG_RESULT) {
154 if (results.size() != 0) {
155 DataLog<dcomplex, dcomplex> logger(getId(), format("stripe[{0}]", stripe - xbegin), "neff", "det");
156 std::string msg = "Found vertical effective indices at: ";
157 for (auto z : results) {
158 msg += str(z) + ", ";
159 logger(z, detS1(z, nrCache[stripe]));
160 }
161 writelog(LOG_RESULT, msg.substr(0, msg.length() - 2));
162 } else
163 writelog(LOG_RESULT, "Did not find any vertical effective indices");
164 }
165
166 return results;
167}
168
169size_t EffectiveIndex2D::findMode(dcomplex neff, Symmetry symmetry) {
170 writelog(LOG_INFO, "Searching for the mode starting from Neff = {0}", str(neff));
171 stageOne();
172 Mode mode(this, symmetry);
173 mode.neff = RootDigger::get(
174 this, [this, &mode](const dcomplex& x) { return this->detS(x, mode); }, log_value, root)
175 ->find(neff);
176 return insertMode(mode);
177}
178
179std::vector<size_t> EffectiveIndex2D::findModes(dcomplex neff1,
180 dcomplex neff2,
181 Symmetry symmetry,
182 size_t resteps,
183 size_t imsteps,
184 dcomplex eps) {
185 stageOne();
186
187 if (eps.imag() == 0.) eps.imag(eps.real());
188
189 if (real(eps) <= 0. || imag(eps) <= 0.) throw BadInput(this->getId(), "bad precision specified");
190
191 double re0 = real(neff1), im0 = imag(neff1);
192 double re1 = real(neff2), im1 = imag(neff2);
193 if (re0 > re1) std::swap(re0, re1);
194 if (im0 > im1) std::swap(im0, im1);
195
196 if (re0 == 0. && re1 == 0.) {
197 re0 = 1e30;
198 re1 = -1e30;
199 for (size_t i = xbegin; i != xend; ++i) {
200 dcomplex n = sqrt(epsilons[i]);
201 if (n.real() < re0) re0 = n.real();
202 if (n.real() > re1) re1 = n.real();
203 }
204 } else if (re0 == 0. || re1 == 0.)
205 throw BadInput(getId(), "bad area to browse specified");
206 if (im0 == 0. && im1 == 0.) {
207 im0 = 1e30;
208 im1 = -1e30;
209 for (size_t i = xbegin; i != xend; ++i) {
210 dcomplex n = sqrt(epsilons[i]);
211 if (n.imag() < im0) im0 = n.imag();
212 if (n.imag() > im1) im1 = n.imag();
213 }
214 }
215 neff1 = dcomplex(re0, im0);
216 neff2 = dcomplex(re1, im1);
217
218 Mode mode(this, symmetry);
219 auto results = findZeros(
220 this, [this, &mode](dcomplex z) { return this->detS(z, mode); }, neff1, neff2, resteps, imsteps, eps);
221
222 std::vector<size_t> idx(results.size());
223
224 if (results.size() != 0) {
225 DataLog<dcomplex, dcomplex> logger(getId(), "Neffs", "Neff", "det");
226 auto refine = RootDigger::get(
227 this, [this, &mode](const dcomplex& v) { return this->detS(v, mode); }, logger, root);
228 std::string msg = "Found modes at: ";
229 for (auto zz : results) {
230 dcomplex z;
231 try {
232 z = refine->find(0.5 * (zz.first + zz.second));
233 } catch (ComputationError&) {
234 continue;
235 }
236 mode.neff = z;
237 idx.push_back(insertMode(mode));
238 msg += str(z) + ", ";
239 }
240 writelog(LOG_RESULT, msg.substr(0, msg.length() - 2));
241 } else
242 writelog(LOG_RESULT, "Did not find any modes");
243
244 return idx;
245}
246
247size_t EffectiveIndex2D::setMode(dcomplex neff, Symmetry sym) {
248 stageOne();
249 Mode mode(this, sym);
250 mode.neff = neff;
251 double det = abs(detS(neff, mode));
252 if (det > root.tolf_max) writelog(LOG_WARNING, "Provided effective index does not correspond to any mode (det = {0})", det);
253 writelog(LOG_INFO, "Setting mode at {0}", str(neff));
254 return insertMode(mode);
255}
256
257void EffectiveIndex2D::onInitialize() {
258 if (!geometry) throw NoGeometryException(getId());
259
260 // Set default mesh
261 if (!mesh) setSimpleMesh();
262
263 xbegin = 0;
264 ybegin = 0;
265 xend = mesh->axis[0]->size() + 1;
266 yend = mesh->axis[1]->size() + 1;
267
268 if (geometry->isExtended(Geometry::DIRECTION_TRAN, false) &&
269 abs(mesh->axis[0]->at(0) - geometry->getChild()->getBoundingBox().lower.c0) < SMALL)
270 xbegin = 1;
271 if (geometry->isExtended(Geometry::DIRECTION_VERT, false) &&
272 abs(mesh->axis[1]->at(0) - geometry->getChild()->getBoundingBox().lower.c1) < SMALL)
273 ybegin = 1;
274 if (geometry->isExtended(Geometry::DIRECTION_TRAN, true) &&
275 abs(mesh->axis[0]->at(mesh->axis[0]->size() - 1) - geometry->getChild()->getBoundingBox().upper.c0) < SMALL)
276 --xend;
277 if (geometry->isExtended(Geometry::DIRECTION_VERT, true) &&
278 abs(mesh->axis[1]->at(mesh->axis[1]->size() - 1) - geometry->getChild()->getBoundingBox().upper.c1) < SMALL)
279 --yend;
280
281 // Assign space for refractive indices cache and stripe effective indices
282 nrCache.assign(xend, std::vector<dcomplex, aligned_allocator<dcomplex>>(yend));
283 epsilons.resize(xend);
284
285 yfields.resize(yend);
286
287 need_gain = false;
288}
289
290void EffectiveIndex2D::onInvalidate() {
291 if (!modes.empty()) {
292 writelog(LOG_DETAIL, "Clearing computed modes");
293 modes.clear();
294 outNeff.fireChanged();
295 outLightMagnitude.fireChanged();
296 outLightE.fireChanged();
297 }
298 recompute_neffs = true;
299}
300
301/********* Here are the computations *********/
302
303void EffectiveIndex2D::updateCache() {
304 bool fresh = initCalculation();
305
306 if (fresh) {
307 // we need to update something
308
309 if (geometry->isSymmetric(Geometry::DIRECTION_TRAN)) {
310 if (fresh) // Make sure we have only positive points
311 for (auto x : *mesh->axis[0])
312 if (x < 0.) throw BadMesh(getId(), "for symmetric geometry no horizontal points can be negative");
313 if (mesh->axis[0]->at(0) == 0.) xbegin = 1;
314 }
315
316 if (!modes.empty()) writelog(LOG_DETAIL, "Clearing computed modes");
317 modes.clear();
318
319 double w = real(2e3 * PI / k0);
320
321 shared_ptr<OrderedAxis> axis0, axis1;
322 {
323 shared_ptr<RectangularMesh<2>> midmesh = mesh->getElementMesh();
324 axis0 = plask::make_shared<OrderedAxis>(*midmesh->axis[0]);
325 axis1 = plask::make_shared<OrderedAxis>(*midmesh->axis[1]);
326 }
327
328 if (xbegin == 0) {
329 if (geometry->isSymmetric(Geometry::DIRECTION_TRAN))
330 axis0->addPoint(0.5 * mesh->axis[0]->at(0));
331 else
332 axis0->addPoint(mesh->axis[0]->at(0) - 2. * OrderedAxis::MIN_DISTANCE);
333 }
334 if (xend == mesh->axis[0]->size() + 1)
335 axis0->addPoint(mesh->axis[0]->at(mesh->axis[0]->size() - 1) + 2. * OrderedAxis::MIN_DISTANCE);
336 if (ybegin == 0) axis1->addPoint(mesh->axis[1]->at(0) - 2. * OrderedAxis::MIN_DISTANCE);
337 if (yend == mesh->axis[1]->size() + 1)
338 axis1->addPoint(mesh->axis[1]->at(mesh->axis[1]->size() - 1) + 2. * OrderedAxis::MIN_DISTANCE);
339
340 writelog(LOG_DEBUG, "Updating refractive indices cache");
341 auto midmesh = plask::make_shared<RectangularMesh<2>>(axis0, axis1, mesh->getIterationOrder());
342 auto temp = inTemperature.hasProvider() ? inTemperature(midmesh) : LazyData<double>(midmesh->size(), 300.);
343 bool have_gain = false;
344 LazyData<Tensor2<double>> gain;
345 auto carriers = inCarriersConcentration.hasProvider() ? inCarriersConcentration(CarriersConcentration::MAJORITY, midmesh)
346 : LazyData<double>(midmesh->size(), 0.);
347
348 for (size_t ix = xbegin; ix < xend; ++ix) {
349 for (size_t iy = ybegin; iy < yend; ++iy) {
350 size_t idx = midmesh->index(ix - xbegin, iy - ybegin);
351 double T = temp[idx];
352 double cc = carriers[idx];
353 auto point = midmesh->at(idx);
354 auto roles = geometry->getRolesAt(point);
355 if (roles.find("QW") == roles.end() && roles.find("QD") == roles.end() && roles.find("gain") == roles.end())
356 nrCache[ix][iy] = geometry->getMaterial(point)->Nr(w, T, cc);
357 else { // we ignore the material absorption as it should be considered in the gain already
358 need_gain = true;
359 if (!have_gain) {
360 gain = inGain(midmesh, w);
361 have_gain = true;
362 }
363 double g = (polarization == TM) ? gain[idx].c11 : gain[idx].c00;
364 nrCache[ix][iy] = dcomplex(real(geometry->getMaterial(point)->Nr(w, T, cc)), w * g * (0.25e-7 / PI));
365 }
366 }
367 }
368 recompute_neffs = true;
369 }
370}
371
372void EffectiveIndex2D::stageOne() {
373 updateCache();
374
375 if (recompute_neffs) {
376 // Compute effective index of the main stripe
377 size_t stripe = mesh->tran()->findIndex(stripex);
378 if (stripe < xbegin)
379 stripe = xbegin;
380 else if (stripe >= xend)
381 stripe = xend - 1;
382 writelog(LOG_DETAIL, "Computing effective index for vertical stripe {0} (polarization {1})", stripe - xbegin,
383 (polarization == TE) ? "TE" : "TM");
384#ifndef NDEBUG
385 {
386 std::stringstream nrs;
387 for (ptrdiff_t j = yend - 1; j >= ptrdiff_t(ybegin); --j) nrs << ", " << str(nrCache[stripe][j]);
388 writelog(LOG_DEBUG, "Nr[{0}] = [{1} ]", stripe - xbegin, nrs.str().substr(1));
389 }
390#endif
391 DataLog<dcomplex, dcomplex> log_stripe(getId(), format("stripe[{0}]", stripe - xbegin), "neff", "det");
392 auto rootdigger = RootDigger::get(
393 this, [&](const dcomplex& x) { return this->detS1(x, nrCache[stripe]); }, log_stripe, stripe_root);
394 if (vneff == 0.) {
395 dcomplex maxn = *std::max_element(nrCache[stripe].begin(), nrCache[stripe].end(),
396 [](const dcomplex& a, const dcomplex& b) { return real(a) < real(b); });
397 vneff = 0.999 * real(maxn);
398 }
399 vneff = rootdigger->find(vneff);
400
401 // Compute field weights
402 computeWeights(stripe);
403
404 // Compute effective indices
405 for (size_t i = xbegin; i < xend; ++i) {
406 epsilons[i] = vneff * vneff;
407 for (size_t j = ybegin; j < yend; ++j) {
408 epsilons[i] += yweights[j] * (nrCache[i][j] * nrCache[i][j] - nrCache[stripe][j] * nrCache[stripe][j]);
409 }
410 }
411 if (maxLoglevel >= LOG_DEBUG) {
412 std::stringstream nrs;
413 for (size_t i = xbegin; i < xend; ++i) {
414 dcomplex n = sqrt(epsilons[i]);
415 if (abs(n.real()) < 1e-10) n.real(0.);
416 if (abs(n.imag()) < 1e-10) n.imag(0.);
417 nrs << ", " << str(n);
418 }
419 writelog(LOG_DEBUG, "vertical neffs = [{0} ]", nrs.str().substr(1));
420 }
421 double rmin = INFINITY, rmax = -INFINITY, imin = INFINITY, imax = -INFINITY;
422 for (size_t i = xbegin; i < xend; ++i) {
423 dcomplex n = sqrt(epsilons[i]);
424 if (real(n) < rmin) rmin = real(n);
425 if (real(n) > rmax) rmax = real(n);
426 if (imag(n) < imin) imin = imag(n);
427 if (imag(n) > imax) imax = imag(n);
428 }
429 writelog(LOG_DETAIL, "Effective index should be between {0} and {1}", str(dcomplex(rmin, imin)), str(dcomplex(rmax, imax)));
430 recompute_neffs = false;
431 }
432}
433
434dcomplex EffectiveIndex2D::detS1(const plask::dcomplex& x,
435 const std::vector<dcomplex, aligned_allocator<dcomplex>>& NR,
436 bool save) {
437 if (save) yfields[ybegin] = Field(0., 1.);
438
439 std::vector<dcomplex, aligned_allocator<dcomplex>> ky(yend);
440 for (size_t i = ybegin; i < yend; ++i) {
441 ky[i] = k0 * sqrt(NR[i] * NR[i] - x * x);
442 if (imag(ky[i]) > 0.) ky[i] = -ky[i];
443 }
444
445 // dcomplex s1 = 1., s2 = 0., s3 = 0., s4 = 1.; // matrix S
446 //
447 // dcomplex phas = 1.;
448 // if (ybegin != 0)
449 // phas = exp(I * ky[ybegin] * (mesh->axis[1]->at(ybegin)-mesh->axis[1]->at(ybegin-1)));
450 //
451 // for (size_t i = ybegin+1; i < yend; ++i) {
452 // // Compute shift inside one layer
453 // s1 *= phas;
454 // s3 *= phas * phas;
455 // s4 *= phas;
456 // // Compute matrix after boundary
457 // dcomplex f = (polarization==TM)? NR[i-1]/NR[i] : 1.;
458 // dcomplex p = 0.5 + 0.5 * ky[i] / ky[i-1] * f*f;
459 // dcomplex m = 1.0 - p;
460 // dcomplex chi = 1. / (p - m * s3);
461 // // F0 = [ (-m*m + p*p)*chi*s1 m*s1*s4*chi + s2 ] [ F2 ]
462 // // B2 = [ (-m + p*s3)*chi s4*chi ] [ B0 ]
463 // s2 += s1*m*chi*s4;
464 // s1 *= (p*p - m*m) * chi;
465 // s3 = (p*s3-m) * chi;
466 // s4 *= chi;
467 // // Compute phase shift for the next step
468 // if (i != mesh->axis[1]->size())
469 // phas = exp(I * ky[i] * (mesh->axis[1]->at(i)-mesh->axis[1]->at(i-1)));
470 //
471 // // Compute fields
472 // if (save) {
473 // dcomplex F = -s2/s1, B = (s1*s4-s2*s3)/s1; // Assume F0 = 0 B0 = 1
474 // double aF = abs(F), aB = abs(B);
475 // // zero very small fields to avoid errors in plotting for long layers
476 // if (aF < 1e-8 * aB) F = 0.;
477 // if (aB < 1e-8 * aF) B = 0.;
478 // yfields[i] = Field(F, B);
479 // }
480 // }
481
482 Matrix T = Matrix::eye();
483 for (size_t i = ybegin; i < yend - 1; ++i) {
484 double d;
485 if (i != ybegin || ybegin != 0)
486 d = mesh->axis[1]->at(i) - mesh->axis[1]->at(i - 1);
487 else
488 d = 0.;
489 dcomplex phas = exp(-I * ky[i] * d);
490 // Transfer through boundary
491 dcomplex f = (polarization == TM) ? (NR[i + 1] / NR[i]) : 1.;
492 dcomplex n = 0.5 * ky[i] / ky[i + 1] * f * f;
493 Matrix T1 = Matrix((0.5 + n), (0.5 - n), (0.5 - n), (0.5 + n));
494 T1.ff *= phas;
495 T1.fb /= phas;
496 T1.bf *= phas;
497 T1.bb /= phas;
498 T = T1 * T;
499 if (save) {
500 dcomplex F = T.fb, B = T.bb; // Assume F0 = 0 B0 = 1
501 double aF = abs(F), aB = abs(B);
502 // zero very small fields to avoid errors in plotting for long layers
503 if (aF < 1e-8 * aB) F = 0.;
504 if (aB < 1e-8 * aF) B = 0.;
505 yfields[i + 1] = Field(F, B);
506 }
507 }
508
509 if (save) {
510 yfields[yend - 1].B = 0.;
511#ifndef NDEBUG
512 std::stringstream nrs;
513 for (size_t i = ybegin; i < yend; ++i) nrs << "), (" << str(yfields[i].F) << ":" << str(yfields[i].B);
514 writelog(LOG_DEBUG, "vertical fields = [{0}) ]", nrs.str().substr(2));
515#endif
516 }
517
518 // return s1*s4 - s2*s3;
519
520 return T.bb; // F0 = 0 Bn = 0
521}
522
523void EffectiveIndex2D::computeWeights(size_t stripe) {
524 // Compute fields
525 detS1(vneff, nrCache[stripe], true);
526
527 yweights.resize(yend);
528 {
529 double ky = abs(imag(k0 * sqrt(nrCache[stripe][ybegin] * nrCache[stripe][ybegin] - vneff * vneff)));
530 yweights[ybegin] = abs2(yfields[ybegin].B) * 0.5 / ky;
531 }
532 {
533 double ky = abs(imag(k0 * sqrt(nrCache[stripe][yend - 1] * nrCache[stripe][yend - 1] - vneff * vneff)));
534 yweights[yend - 1] = abs2(yfields[yend - 1].F) * 0.5 / ky;
535 }
536 double sum = yweights[ybegin] + yweights[yend - 1];
537
538 for (size_t i = ybegin + 1; i < yend - 1; ++i) {
539 double d = mesh->axis[1]->at(i) - mesh->axis[1]->at(i - 1);
540 dcomplex ky = k0 * sqrt(nrCache[stripe][i] * nrCache[stripe][i] - vneff * vneff);
541 if (imag(ky) > 0.) ky = -ky;
542 dcomplex w_ff, w_bb, w_fb, w_bf;
543 if (d != 0.) {
544 if (abs(imag(ky)) > SMALL) {
545 dcomplex kk = ky - conj(ky);
546 w_ff = (exp(-I * d * kk) - 1.) / kk;
547 w_bb = -(exp(+I * d * kk) - 1.) / kk;
548 } else
549 w_ff = w_bb = dcomplex(0., -d);
550 if (abs(real(ky)) > SMALL) {
551 dcomplex kk = ky + conj(ky);
552 w_fb = (exp(-I * d * kk) - 1.) / kk;
553 w_bf = -(exp(+I * d * kk) - 1.) / kk;
554 } else
555 w_ff = w_bb = dcomplex(0., -d);
556 dcomplex weight = yfields[i].F * conj(yfields[i].F) * w_ff + yfields[i].F * conj(yfields[i].B) * w_fb +
557 yfields[i].B * conj(yfields[i].F) * w_bf + yfields[i].B * conj(yfields[i].B) * w_bb;
558 yweights[i] = -imag(weight);
559 } else
560 yweights[i] = 0.;
561 sum += yweights[i];
562 }
563
564 sum = 1. / sum;
565 double fact = sqrt(sum);
566 for (size_t i = ybegin; i < yend; ++i) {
567 yweights[i] *= sum;
568 yfields[i] *= fact;
569 }
570 // #ifndef NDEBUG
571 // {
572 // std::stringstream nrs; for (size_t i = ybegin; i < yend; ++i) nrs << ", " << str(weights[i]);
573 // writelog(LOG_DEBUG, "vertical weights = [{0}) ]", nrs.str().substr(2));
574 // }
575 // #endif
576}
577
578void EffectiveIndex2D::normalizeFields(Mode& mode, const std::vector<dcomplex, aligned_allocator<dcomplex>>& kx) {
579 size_t start;
580
581 double sum = mode.xweights[xend - 1] = abs2(mode.xfields[xend - 1].F) * 0.5 / abs(imag(kx[xend - 1]));
582 if (mode.symmetry == SYMMETRY_NONE) {
583 sum += mode.xweights[0] = abs2(mode.xfields[xbegin].B) * 0.5 / abs(imag(kx[xbegin]));
584 start = xbegin + 1;
585 } else {
586 start = xbegin;
587 }
588
589 for (size_t i = start; i < xend - 1; ++i) {
590 double d = mesh->axis[0]->at(i) - ((i == 0) ? 0. : mesh->axis[0]->at(i - 1));
591 dcomplex w_ff, w_bb, w_fb, w_bf;
592 if (d != 0.) {
593 if (abs(imag(kx[i])) > SMALL) {
594 dcomplex kk = kx[i] - conj(kx[i]);
595 w_ff = (exp(-I * d * kk) - 1.) / kk;
596 w_bb = -(exp(+I * d * kk) - 1.) / kk;
597 } else
598 w_ff = w_bb = dcomplex(0., -d);
599 if (abs(real(kx[i])) > SMALL) {
600 dcomplex kk = kx[i] + conj(kx[i]);
601 w_fb = (exp(-I * d * kk) - 1.) / kk;
602 w_bf = -(exp(+I * d * kk) - 1.) / kk;
603 } else
604 w_ff = w_bb = dcomplex(0., -d);
605 mode.xweights[i] =
606 -imag(mode.xfields[i].F * conj(mode.xfields[i].F) * w_ff + mode.xfields[i].F * conj(mode.xfields[i].B) * w_fb +
607 mode.xfields[i].B * conj(mode.xfields[i].F) * w_bf + mode.xfields[i].B * conj(mode.xfields[i].B) * w_bb);
608 sum += mode.xweights[i];
609 }
610 }
611 if (mode.symmetry != SYMMETRY_NONE) sum *= 2.;
612
613 // Consider loss on the mirror
614 double R1, R2;
615 if (mirrors) {
616 std::tie(R1, R2) = *mirrors;
617 } else {
618 const double lambda = real(2e3 * PI / k0);
619 const double n = real(mode.neff);
620 const double n1 = real(geometry->getFrontMaterial()->Nr(lambda, 300.)),
621 n2 = real(geometry->getBackMaterial()->Nr(lambda, 300.));
622 R1 = (n - n1) / (n + n1);
623 R1 *= R1;
624 R2 = (n - n2) / (n + n2);
625 R2 *= R2;
626 }
627
628 if (emission == FRONT) {
629 if (R1 == 1.) this->writelog(LOG_WARNING, "Mirror reflection on emission side equal to 1. Field will be infinite.");
630 sum *= (1. - R1);
631 } else {
632 if (R2 == 1.) this->writelog(LOG_WARNING, "Mirror reflection on emission side equal to 1. Field will be infinite.");
633 sum *= (1. - R2);
634 }
635
636 double ff = 1e12 / sum; // 1e12 because intensity in W/m² and integral computed in µm
637 double f = sqrt(ff);
638
639 for (auto& val : mode.xfields) val *= f;
640 for (auto& val : mode.xweights) val *= ff;
641}
642
643double EffectiveIndex2D::getTotalAbsorption(Mode& mode) {
644 if (!mode.have_fields) detS(mode.neff, mode, true);
645
646 double result = 0.;
647
648 for (size_t ix = 0; ix < xend; ++ix) {
649 for (size_t iy = ybegin; iy < yend; ++iy) {
650 double absp = -2. * real(nrCache[ix][iy]) * imag(nrCache[ix][iy]);
651 result += absp * mode.xweights[ix] * yweights[iy]; // [dV] = µm³
652 }
653 }
654 if (mode.symmetry != SYMMETRY_NONE) result *= 2.;
655 result *= 1e-9 * real(k0) * mode.power; // 1e-9: µm³ / nm -> m², ½ is already hidden in mode.power
656 return result;
657}
658
659double EffectiveIndex2D::getTotalAbsorption(std::size_t num) {
660 if (modes.size() <= num) throw NoValue("absorption");
661
662 return getTotalAbsorption(modes[num]);
663}
664
665dcomplex EffectiveIndex2D::detS(const dcomplex& x, EffectiveIndex2D::Mode& mode, bool save) {
666 // Adjust for mirror losses
667 dcomplex neff2 = dcomplex(real(x), imag(x) - getMirrorLosses(x));
668 neff2 *= neff2;
669
670 std::vector<dcomplex, aligned_allocator<dcomplex>> kx(xend);
671 for (size_t i = xbegin; i < xend; ++i) {
672 kx[i] = k0 * sqrt(epsilons[i] - neff2);
673 if (imag(kx[i]) > 0.) kx[i] = -kx[i];
674 }
675
676 aligned_unique_ptr<Matrix[]> matrices;
677 if (save) matrices.reset(aligned_malloc<Matrix>(xend - 1));
678
679 Matrix T = Matrix::eye();
680 for (size_t i = xbegin; i < xend - 1; ++i) {
681 double d;
682 if (i != xbegin)
683 d = mesh->axis[0]->at(i) - mesh->axis[0]->at(i - 1);
684 else if (mode.symmetry != SYMMETRY_NONE)
685 d = mesh->axis[0]->at(i); // we have symmetry, so beginning of the transfer matrix is at the axis
686 else
687 d = 0.;
688 dcomplex phas = exp(-I * kx[i] * d);
689 // Transfer through boundary
690 dcomplex f = (polarization == TE) ? (epsilons[i + 1] / epsilons[i]) : 1.;
691 dcomplex n = 0.5 * kx[i] / kx[i + 1] * f;
692 Matrix T1 = Matrix((0.5 + n), (0.5 - n), (0.5 - n), (0.5 + n));
693 T1.ff *= phas;
694 T1.fb /= phas;
695 T1.bf *= phas;
696 T1.bb /= phas;
697 if (save) matrices[i] = T1;
698 T = T1 * T;
699 }
700
701 if (save) {
702 mode.neff = x;
703
704 mode.xfields[xend - 1] = Field(1., 0.);
705 for (size_t i = xend - 1; i != xbegin; --i) {
706 mode.xfields[i - 1] = matrices[i - 1].solve(mode.xfields[i]);
707 }
708 normalizeFields(mode, kx);
709 mode.have_fields = true;
710#ifndef NDEBUG
711 {
712 std::stringstream nrs;
713 for (size_t i = xbegin; i < xend; ++i) nrs << "), (" << str(mode.xfields[i].F) << ":" << str(mode.xfields[i].B);
714 writelog(LOG_DEBUG, "horizontal fields = [{0}) ]", nrs.str().substr(2));
715 }
716#endif
717 }
718
719 if (mode.symmetry == SYMMETRY_POSITIVE)
720 return T.bf + T.bb; // B0 = F0 Bn = 0
721 else if (mode.symmetry == SYMMETRY_NEGATIVE)
722 return T.bf - T.bb; // B0 = -F0 Bn = 0
723 else
724 return T.bb; // F0 = 0 Bn = 0
725}
726
727template <typename FieldT> struct EffectiveIndex2D::FieldDataBase : public LazyDataImpl<FieldT> {
728 EffectiveIndex2D* solver;
729 std::size_t num;
730 std::vector<dcomplex, aligned_allocator<dcomplex>> kx, ky;
731 std::size_t stripe;
732
733 FieldDataBase(EffectiveIndex2D* solver, std::size_t num)
734 : solver(solver), num(num), kx(solver->xend), ky(solver->yend), stripe(solver->mesh->tran()->findIndex(solver->stripex)) {
735 dcomplex neff = solver->modes[num].neff;
736 if (!solver->modes[num].have_fields) solver->detS(neff, solver->modes[num], true);
737
738 if (stripe < solver->xbegin)
739 stripe = solver->xbegin;
740 else if (stripe >= solver->xend)
741 stripe = solver->xend - 1;
742
743 solver->writelog(LOG_INFO, "Computing field distribution for Neff = {0}", str(neff));
744
745 dcomplex neff2 = dcomplex(real(neff), imag(neff) - solver->getMirrorLosses(real(neff)));
746 neff2 *= neff2;
747 for (size_t i = 0; i < solver->xend; ++i) {
748 kx[i] = solver->k0 * sqrt(solver->epsilons[i] - neff2);
749 if (imag(kx[i]) > 0.) kx[i] = -kx[i];
750 }
751
752 for (size_t i = solver->ybegin; i < solver->yend; ++i) {
753 ky[i] = solver->k0 * sqrt(solver->nrCache[stripe][i] * solver->nrCache[stripe][i] - solver->vneff * solver->vneff);
754 if (imag(ky[i]) > 0.) ky[i] = -ky[i];
755 }
756
757 setScale();
758 }
759
760 protected:
761 inline FieldT value(dcomplex val) const;
762 double scale;
763 void setScale();
764};
765
766template <> void EffectiveIndex2D::FieldDataBase<double>::setScale() { scale = 1e-3 * solver->modes[num].power; }
767
768template <> double EffectiveIndex2D::FieldDataBase<double>::value(dcomplex val) const { return scale * abs2(val); }
769
770template <> void EffectiveIndex2D::FieldDataBase<Vec<3, dcomplex>>::setScale() {
771 // <M> = ½ E conj(E) / Z0
772 scale = sqrt(2e-3 * solver->modes[num].power * phys::Z0);
773}
774
775template <> Vec<3, dcomplex> EffectiveIndex2D::FieldDataBase<Vec<3, dcomplex>>::value(dcomplex val) const {
776 if (solver->getPolarization() == TE)
777 return Vec<3, dcomplex>(0., scale * val, 0.);
778 else
779 return Vec<3, dcomplex>(0., 0., scale * val);
780}
781
782template <typename FieldT> struct EffectiveIndex2D::FieldDataInefficient : public EffectiveIndex2D::FieldDataBase<FieldT> {
783 shared_ptr<const MeshD<2>> dst_mesh;
784
785 FieldDataInefficient(EffectiveIndex2D* solver, std::size_t num, const shared_ptr<const MeshD<2>>& dst_mesh)
786 : FieldDataBase<FieldT>(solver, num), dst_mesh(dst_mesh) {}
787
788 size_t size() const override { return dst_mesh->size(); }
789
790 FieldT at(size_t idx) const override {
791 auto point = dst_mesh->at(idx);
792 double x = point.tran();
793 double y = point.vert();
794
795 bool negate = false;
796 if (x < 0. && this->solver->modes[this->num].symmetry != EffectiveIndex2D::SYMMETRY_NONE) {
797 x = -x;
798 if (this->solver->modes[this->num].symmetry == EffectiveIndex2D::SYMMETRY_NEGATIVE) negate = true;
799 }
800 size_t ix = this->solver->mesh->tran()->findIndex(x);
801 if (ix >= this->solver->xend) ix = this->solver->xend - 1;
802 if (ix < this->solver->xbegin) ix = this->solver->xbegin;
803 if (ix != 0)
804 x -= this->solver->mesh->tran()->at(ix - 1);
805 else if (this->solver->modes[this->num].symmetry == EffectiveIndex2D::SYMMETRY_NONE)
806 x -= this->solver->mesh->tran()->at(0);
807 dcomplex phasx = exp(-I * this->kx[ix] * x);
808 dcomplex val = this->solver->modes[this->num].xfields[ix].F * phasx + this->solver->modes[this->num].xfields[ix].B / phasx;
809 if (negate) val = -val;
810
811 size_t iy = this->solver->mesh->vert()->findIndex(y);
812 if (iy >= this->solver->yend) iy = this->solver->yend - 1;
813 if (iy < this->solver->ybegin) iy = this->solver->ybegin;
814 y -= this->solver->mesh->vert()->at(max(int(iy) - 1, 0));
815 dcomplex phasy = exp(-I * this->ky[iy] * y);
816 val *= this->solver->yfields[iy].F * phasy + this->solver->yfields[iy].B / phasy;
817
818 return this->value(val);
819 }
820};
821
822template <typename FieldT> struct EffectiveIndex2D::FieldDataEfficient : public EffectiveIndex2D::FieldDataBase<FieldT> {
823 shared_ptr<const RectangularMesh<2>> rect_mesh;
824 std::vector<dcomplex, aligned_allocator<dcomplex>> valx, valy;
825
826 size_t size() const override { return rect_mesh->size(); }
827
828 FieldDataEfficient(EffectiveIndex2D* solver, std::size_t num, const shared_ptr<const RectangularMesh<2>>& rect_mesh)
829 : FieldDataBase<FieldT>(solver, num),
830 rect_mesh(rect_mesh),
831 valx(rect_mesh->tran()->size()),
832 valy(rect_mesh->vert()->size()) {
833PLASK_OMP_PARALLEL
834 {
835#pragma omp for nowait
836 for (plask::openmp_size_t idx = 0; idx < rect_mesh->tran()->size(); ++idx) {
837 double x = rect_mesh->tran()->at(idx);
838 bool negate = false;
839 if (x < 0. && this->solver->modes[num].symmetry != EffectiveIndex2D::SYMMETRY_NONE) {
840 x = -x;
841 if (this->solver->modes[num].symmetry == EffectiveIndex2D::SYMMETRY_NEGATIVE) negate = true;
842 }
843 size_t ix = this->solver->mesh->tran()->findIndex(x);
844 if (ix >= this->solver->xend) ix = this->solver->xend - 1;
845 if (ix < this->solver->xbegin) ix = this->solver->xbegin;
846 if (ix != 0)
847 x -= this->solver->mesh->tran()->at(ix - 1);
848 else if (this->solver->modes[num].symmetry == EffectiveIndex2D::SYMMETRY_NONE)
849 x -= this->solver->mesh->tran()->at(0);
850 dcomplex phasx = exp(-I * this->kx[ix] * x);
851 dcomplex val =
852 this->solver->modes[this->num].xfields[ix].F * phasx + this->solver->modes[num].xfields[ix].B / phasx;
853 if (negate) val = -val;
854 valx[idx] = val;
855 }
856
857#pragma omp for
858 for (plask::openmp_size_t idy = 0; idy < rect_mesh->vert()->size(); ++idy) {
859 double y = rect_mesh->vert()->at(idy);
860 size_t iy = this->solver->mesh->vert()->findIndex(y);
861 if (iy >= this->solver->yend) iy = this->solver->yend - 1;
862 if (iy < this->solver->ybegin) iy = this->solver->ybegin;
863 y -= solver->mesh->vert()->at(max(int(iy) - 1, 0));
864 dcomplex phasy = exp(-I * this->ky[iy] * y);
865 valy[idy] = this->solver->yfields[iy].F * phasy + this->solver->yfields[iy].B / phasy;
866 }
867 }
868 // Free no longer needed memory
869 this->kx.clear();
870 this->ky.clear();
871 }
872
873 FieldT at(size_t idx) const override {
874 size_t i0 = rect_mesh->index0(idx);
875 size_t i1 = rect_mesh->index1(idx);
876 return this->value(valx[i0] * valy[i1]);
877 }
878
879 DataVector<const FieldT> getAll() const override {
880 DataVector<FieldT> results(rect_mesh->size());
881 if (rect_mesh->getIterationOrder() == RectangularMesh<2>::ORDER_10) {
882PLASK_OMP_PARALLEL_FOR
883 for (plask::openmp_size_t i1 = 0; i1 < rect_mesh->axis[1]->size(); ++i1) {
884 FieldT* data = results.data() + i1 * rect_mesh->axis[0]->size();
885 for (size_t i0 = 0; i0 < rect_mesh->axis[0]->size(); ++i0) {
886 dcomplex f = valx[i0] * valy[i1];
887 data[i0] = this->value(f);
888 }
889 }
890 } else {
891PLASK_OMP_PARALLEL_FOR
892 for (plask::openmp_size_t i0 = 0; i0 < rect_mesh->axis[0]->size(); ++i0) {
893 FieldT* data = results.data() + i0 * rect_mesh->axis[1]->size();
894 for (size_t i1 = 0; i1 < rect_mesh->axis[1]->size(); ++i1) {
895 dcomplex f = valx[i0] * valy[i1];
896 data[i1] = this->value(f);
897 }
898 }
899 }
900 return results;
901 }
902};
903
904const LazyData<double> EffectiveIndex2D::getLightMagnitude(std::size_t num,
905 shared_ptr<const plask::MeshD<2>> dst_mesh,
906 plask::InterpolationMethod) {
907 this->writelog(LOG_DEBUG, "Getting light intensity");
908
909 if (auto rect_mesh = dynamic_pointer_cast<const RectangularMesh<2>>(dst_mesh))
910 return LazyData<double>(new FieldDataEfficient<double>(this, num, rect_mesh));
911 else
912 return LazyData<double>(new FieldDataInefficient<double>(this, num, dst_mesh));
913}
914
915const LazyData<Vec<3, dcomplex>> EffectiveIndex2D::getElectricField(std::size_t num,
916 shared_ptr<const plask::MeshD<2>> dst_mesh,
917 plask::InterpolationMethod) {
918 this->writelog(LOG_DEBUG, "Getting optical electric field");
919
920 if (auto rect_mesh = dynamic_pointer_cast<const RectangularMesh<2>>(dst_mesh))
921 return LazyData<Vec<3, dcomplex>>(new FieldDataEfficient<Vec<3, dcomplex>>(this, num, rect_mesh));
922 else
923 return LazyData<Vec<3, dcomplex>>(new FieldDataInefficient<Vec<3, dcomplex>>(this, num, dst_mesh));
924}
925
926const LazyData<dcomplex> EffectiveIndex2D::getRefractiveIndex(RefractiveIndex::EnumType,
927 shared_ptr<const MeshD<2>> dst_mesh,
928 dcomplex lam,
929 InterpolationMethod) {
930 if (!isnan(real(lam))) throw BadInput(this->getId(), "wavelength cannot be specified for outRefractiveIndex in this solver");
931 this->writelog(LOG_DEBUG, "Getting refractive indices");
932 updateCache();
933 InterpolationFlags flags(geometry);
934 return LazyData<dcomplex>(dst_mesh->size(), [this, dst_mesh, flags](size_t i) -> dcomplex {
935 auto point = flags.wrap(dst_mesh->at(i));
936 size_t ix = this->mesh->axis[0]->findIndex(point.c0);
937 if (ix < this->xbegin) ix = this->xbegin;
938 size_t iy = this->mesh->axis[1]->findIndex(point.c1);
939 return this->nrCache[ix][iy];
940 });
941}
942
943struct EffectiveIndex2D::HeatDataImpl : public LazyDataImpl<double> {
944 EffectiveIndex2D* solver;
945 shared_ptr<const MeshD<2>> dest_mesh;
946 InterpolationFlags flags;
947 std::vector<LazyData<double>> EE;
948 dcomplex lam0;
949
950 HeatDataImpl(EffectiveIndex2D* solver, const shared_ptr<const MeshD<2>>& dst_mesh, InterpolationMethod method)
951 : solver(solver), dest_mesh(dst_mesh), flags(solver->geometry), EE(solver->modes.size()), lam0(2e3 * PI / solver->k0) {
952 for (std::size_t m = 0; m != solver->modes.size(); ++m) EE[m] = solver->getLightMagnitude(m, dst_mesh, method);
953 }
954
955 size_t size() const override { return dest_mesh->size(); }
956
957 double at(size_t j) const override {
958 double result = 0.;
959 auto point = flags.wrap(dest_mesh->at(j));
960 size_t ix = solver->mesh->axis[0]->findIndex(point.c0);
961 if (ix < solver->xbegin) ix = solver->xbegin;
962 size_t iy = solver->mesh->axis[1]->findIndex(point.c1);
963 for (size_t m = 0; m != solver->modes.size(); ++m) { // we sum heats from all modes
964 result += EE[m][j]; // 1e9: 1/nm -> 1/m
965 }
966 double absp = -2. * real(solver->nrCache[ix][iy]) * imag(solver->nrCache[ix][iy]);
967 result *= 1e6 * real(solver->k0) * absp;
968 return result;
969 }
970};
971
972const LazyData<double> EffectiveIndex2D::getHeat(shared_ptr<const MeshD<2>> dst_mesh, plask::InterpolationMethod method) {
973 // This is somehow naive implementation using the field value from the mesh points. The heat may be slightly off
974 // in case of fast varying light intensity and too sparse mesh.
975 writelog(LOG_DEBUG, "Getting heat absorbed from {0} mode{1}", modes.size(), (modes.size() == 1) ? "" : "s");
976 if (modes.size() == 0) return LazyData<double>(dst_mesh->size(), 0.);
977 return LazyData<double>(new HeatDataImpl(this, dst_mesh, method));
978}
979
980}}} // namespace plask::optical::effective