PLaSK library

Loading...
Searching...
No Matches
ferminew.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 "ferminew.hpp"
15
16 namespace plask { namespace solvers { namespace FermiNew {
17
18 using phys::nm_to_eV;
19
20 template <typename T> struct ptrVector {
21 std::vector<T*> data;
22 ~ptrVector() {
23 for (T* ptr : data) delete ptr;
24 }
25};
26
27 template <typename GeometryType>
28 FermiNewGainSolver<GeometryType>::FermiNewGainSolver(const std::string& name)
29 : SolverWithMesh<GeometryType, MeshAxis>(name),
30 outGain(this, &FermiNewGainSolver<GeometryType>::getGain),
31 outLuminescence(this, &FermiNewGainSolver<GeometryType>::getLuminescence) {
32 Tref = NAN; // (K), only for this temperature energy levels are calculated
33 inTemperature = 300.; // temperature receiver has some sensible value
34 condQWshift = 0.; // (eV)
35 valeQWshift = 0.; // (eV)
36 QWwidthMod = 100.; // (-) (if equal to 10 - differences in QW widths are to big)
37 roughness = 1.00; // (-)
38 lifetime = 0.0; // [ps]
39 matrixElem = 0.; // [m0*eV]
40 differenceQuotient = 0.01; // [%]
41 strains = false;
42 build_struct_once = true;
43 adjust_widths = true;
44 inTemperature.changedConnectMethod(this, &FermiNewGainSolver<GeometryType>::onInputChange);
45 inCarriersConcentration.changedConnectMethod(this, &FermiNewGainSolver<GeometryType>::onInputChange);
46}
47
48 template <typename GeometryType> FermiNewGainSolver<GeometryType>::~FermiNewGainSolver() {
49 disconnectModGeometry();
50 inTemperature.changedDisconnectMethod(this, &FermiNewGainSolver<GeometryType>::onInputChange);
51 inCarriersConcentration.changedDisconnectMethod(this, &FermiNewGainSolver<GeometryType>::onInputChange);
52}
53
54 template <typename GeometryType> void FermiNewGainSolver<GeometryType>::loadConfiguration(XMLReader& reader, Manager& manager) {
55 // Load a configuration parameter from XML.
56 while (reader.requireTagOrEnd()) {
57 std::string param = reader.getNodeName();
58 if (param == "config") {
59 roughness = reader.getAttribute<double>("roughness", roughness);
60 lifetime = reader.getAttribute<double>("lifetime", lifetime);
61 matrixElem = reader.getAttribute<double>("matrix-elem", matrixElem);
62 condQWshift = reader.getAttribute<double>("cond-qw-shift", condQWshift);
63 valeQWshift = reader.getAttribute<double>("vale-qw-shift", valeQWshift);
64 Tref = reader.getAttribute<double>("Tref", Tref);
65 strains = reader.getAttribute<bool>("strained", strains);
66 adjust_widths = reader.getAttribute<bool>("adjust-layers", adjust_widths);
67 build_struct_once = reader.getAttribute<bool>("fast-levels", build_struct_once);
68 if (reader.hasAttribute("substrate")) {
69 substrateMaterial = MaterialsDB::getDefault().get(reader.requireAttribute<std::string>("substrate"));
70 explicitSubstrate = true;
71 }
72 reader.requireTagEnd();
73 // } else if (param == "levels") {
74 // std::string els, hhs, lhs;
75 // if (reader.hasAttribute("el") || reader.hasAttribute("hh") || reader.hasAttribute("lh")) {
76 // els = reader.requireAttribute("el");
77 // hhs = reader.requireAttribute("hh");
78 // lhs = reader.requireAttribute("lh");
79 // reader.requireTagEnd();
80 // } else {
81 // while (reader.requireTagOrEnd()) {
82 // if (reader.getNodeName() == "el") els = reader.requireTextInCurrentTag();
83 // else if (reader.getNodeName() == "hh") hhs = reader.requireTextInCurrentTag();
84 // else if (reader.getNodeName() == "lh") lhs = reader.requireTextInCurrentTag();
85 // else throw XMLUnexpectedElementException(reader, "<el>, <hh>, or <lh>");
86 // }
87 // if (els == "") throw XMLUnexpectedElementException(reader, "<el>");
88 // if (hhs == "") throw XMLUnexpectedElementException(reader, "<hh>");
89 // if (lhs == "") throw XMLUnexpectedElementException(reader, "<lh>");
90 // }
91 // boost::char_separator<char> sep(", ");
92 // boost::tokenizer<boost::char_separator<char>> elt(els, sep), hht(hhs, sep), lht(lhs, sep);
93 // /*double *el = nullptr, *hh = nullptr, *lh = nullptr;
94 // try {
95 // el = new double[std::distance(elt.begin(), elt.end())+1];
96 // hh = new double[std::distance(hht.begin(), hht.end())+1];
97 // lh = new double[std::distance(lht.begin(), lht.end())+1];
98 // double* e = el; for (const auto& i: elt) *(e++) = - boost::lexical_cast<double>(i); *e = 1.;
99 // double* h = hh; for (const auto& i: hht) *(h++) = - boost::lexical_cast<double>(i); *h = 1.;
100 // double* l = lh; for (const auto& i: lht) *(l++) = - boost::lexical_cast<double>(i); *l = 1.;
101 // } catch(...) {
102 // delete[] el; delete[] hh; delete[] lh;
103 // }*/
104 // std::unique_ptr<double[]> el(new double[std::distance(elt.begin(), elt.end())+1]);
105 // std::unique_ptr<double[]> hh(new double[std::distance(hht.begin(), hht.end())+1]);
106 // std::unique_ptr<double[]> lh(new double[std::distance(lht.begin(), lht.end())+1]);
107 // double* e = el.get(); for (const auto& i: elt) *(e++) = - boost::lexical_cast<double>(i); *e = 1.;
108 // double* h = hh.get(); for (const auto& i: hht) *(h++) = - boost::lexical_cast<double>(i); *h = 1.;
109 // double* l = lh.get(); for (const auto& i: lht) *(l++) = - boost::lexical_cast<double>(i); *l = 1.;
110 // /*if (extern_levels) {
111 // delete[] extern_levels->el; delete[] extern_levels->hh; delete[] extern_levels->lh;
112 // }
113 // extern_levels.reset(ExternalLevels(el.release(), hh.release(), lh.release()));*/
114 } else {
115 if (param == "geometry") {
116 auto name = reader.getAttribute("mod");
117 if (name) {
118 auto found = manager.geometrics.find(*name);
119 if (found == manager.geometrics.end())
120 throw BadInput(this->getId(), "geometry '{0}' not found", *name);
121 else {
122 auto geometry = dynamic_pointer_cast<GeometryType>(found->second);
123 if (!geometry) throw BadInput(this->getId(), "geometry '{0}' of wrong type", *name);
124 this->setModGeometry(geometry);
125 }
126 }
127 }
128 this->parseStandardConfiguration(reader, manager, "<geometry>, <mesh>, <levels>, or <config>");
129 }
130 }
131}
132
133 template <typename GeometryType>
134 void FermiNewGainSolver<GeometryType>::onInitialize() // In this function check if geometry and mesh are set
135{
136 if (!this->geometry) throw NoGeometryException(this->getId());
137
138 prepareActiveRegionsInfo();
139 if (build_struct_once) {
140 region_levels.resize(regions.size());
141 }
142
143 outGain.fireChanged();
144 outLuminescence.fireChanged();
145}
146
147 template <typename GeometryType>
148 void FermiNewGainSolver<GeometryType>::onInvalidate() // This will be called when e.g. geometry or mesh changes and
149 // your results become outdated
150{
151 region_levels.clear();
152}
153
154/*template <typename GeometryType>
155 void FermiNewGainSolver<GeometryType>::compute() {
156 this->initCalculation(); // This must be called before any calculation!
157}*/
158
159 template <typename GeometryType>
160std::list<typename FermiNewGainSolver<GeometryType>::ActiveRegionData> FermiNewGainSolver<GeometryType>::detectActiveRegions(
161 const shared_ptr<GeometryType>& geometry) {
162 std::list<typename FermiNewGainSolver<GeometryType>::ActiveRegionData> regions;
163
164 shared_ptr<RectangularMesh<2>> mesh = makeGeometryGrid(geometry->getChild());
165 shared_ptr<RectangularMesh<2>> points = mesh->getElementMesh();
166
167 size_t ileft = 0, iright = points->axis[0]->size();
168 bool in_active = false;
169
170 bool added_bottom_cladding = false;
171 bool added_top_cladding = false;
172
173 for (size_t r = 0; r < points->axis[1]->size(); ++r) {
174 bool had_active = false; // indicates if we had active region in this layer
175 shared_ptr<Material> layer_material;
176 bool layer_QW = false;
177
178 for (size_t c = 0; c < points->axis[0]->size(); ++c) { // In the (possible) active region
179 auto point = points->at(c, r);
180 auto tags = geometry->getRolesAt(point);
181 bool active = false;
182 for (const auto& tag : tags)
183 if (tag.substr(0, 6) == "active") {
184 active = true;
185 break;
186 }
187 bool QW = tags.find("QW") != tags.end() /* || tags.find("QD") != tags.end()*/;
188 bool substrate = tags.find("substrate") != tags.end();
189
190 if (substrate) {
191 if (this->explicitSubstrate)
192 this->writelog(LOG_WARNING, "Explicit substrate layer specified, role 'substrate' ignored");
193 else {
194 if (!substrateMaterial)
195 substrateMaterial = geometry->getMaterial(point);
196 else if (*substrateMaterial != *geometry->getMaterial(point))
197 throw Exception("{0}: Non-uniform substrate layer.", this->getId());
198 }
199 }
200
201 if (QW && !active) throw Exception("{0}: All marked quantum wells must belong to marked active region.", this->getId());
202
203 if (c < ileft) {
204 if (active) throw Exception("{0}: Left edge of the active region not aligned.", this->getId());
205 } else if (c >= iright) {
206 if (active) throw Exception("{0}: Right edge of the active region not aligned.", this->getId());
207 } else {
208 // Here we are inside potential active region
209 if (active) {
210 if (!had_active) {
211 if (!in_active) { // active region is starting set-up new region info
212 regions.emplace_back(mesh->at(c, r));
213 added_bottom_cladding = (tags.find("cladding") != tags.end());
214 added_top_cladding = false;
215 ileft = c;
216 } else {
217 if (added_top_cladding)
218 throw Exception(
219 "{0}: Only the first or the last layer in an active region can have 'cladding' "
220 "role",
221 this->getId());
222 if (tags.find("cladding") != tags.end()) added_top_cladding = true;
223 }
224 layer_material = geometry->getMaterial(point);
225 layer_QW = QW;
226 } else {
227 if (*layer_material != *geometry->getMaterial(point))
228 throw Exception("{0}: Non-uniform active region layer.", this->getId());
229 if (layer_QW != QW)
230 throw Exception("{0}: Quantum-well role of the active region layer not consistent.", this->getId());
231 }
232 } else if (had_active) {
233 if (!in_active) {
234 iright = c;
235 } else
236 throw Exception("{0}: Right edge of the active region not aligned.", this->getId());
237 }
238 had_active |= active;
239 }
240 }
241 in_active = had_active;
242
243 // Now fill-in the layer info
244 if (!regions.empty()) {
245 ActiveRegionData& region = regions.back();
246 if (!added_bottom_cladding) {
247 if (r == 0) throw Exception("{0}: Active region cannot start from the edge of the structure.", this->getId());
248 // add layer below active region (cladding) LUKASZ
249 auto bottom_material = geometry->getMaterial(points->at(ileft, r - 1));
250 for (size_t cc = ileft; cc < iright; ++cc)
251 if (*geometry->getMaterial(points->at(cc, r - 1)) != *bottom_material)
252 throw Exception("{0}: Material below active region not uniform.", this->getId());
253 double w = mesh->axis[0]->at(iright) - mesh->axis[0]->at(ileft);
254 double h = mesh->axis[1]->at(r) - mesh->axis[1]->at(r - 1);
255 region.origin += Vec<2>(0., -h);
256 // this->writelog(LOG_DEBUG, "Adding bottom cladding; h = {0}",h);
257 region.layers->push_back(plask::make_shared<Block<2>>(Vec<2>(w, h), bottom_material));
258 added_bottom_cladding = true;
259 }
260
261 double h = mesh->axis[1]->at(r + 1) - mesh->axis[1]->at(r);
262 double w = mesh->axis[0]->at(iright) - mesh->axis[0]->at(ileft);
263 if (in_active) {
264 size_t n = region.layers->getChildrenCount();
265 shared_ptr<Block<2>> last;
266 if (n > 0)
267 last = static_pointer_cast<Block<2>>(
268 static_pointer_cast<Translation<2>>(region.layers->getChildNo(n - 1))->getChild());
269 assert(!last || last->size.c0 == w);
270 if (last && layer_material == last->getRepresentativeMaterial() && layer_QW == region.isQW(region.size() - 1)) {
271 // TODO check if usage of getRepresentativeMaterial is fine here (was material)
272 last->setSize(w, last->size.c1 + h);
273 } else {
274 auto layer = plask::make_shared<Block<2>>(Vec<2>(w, h), layer_material);
275 if (layer_QW) layer->addRole("QW");
276 region.layers->push_back(layer);
277 // if (layer_QW) this->writelog(LOG_DEBUG, "Adding qw; h = {0}",h);
278 // else this->writelog(LOG_DEBUG, "Adding barrier; h = {0}",h);
279 }
280 } else {
281 if (!added_top_cladding) {
282 // add layer above active region (top cladding)
283 auto top_material = geometry->getMaterial(points->at(ileft, r));
284 for (size_t cc = ileft; cc < iright; ++cc)
285 if (*geometry->getMaterial(points->at(cc, r)) != *top_material)
286 throw Exception("{0}: Material above quantum well not uniform.", this->getId());
287 region.layers->push_back(plask::make_shared<Block<2>>(Vec<2>(w, h), top_material));
288 // this->writelog(LOG_DEBUG, "Adding top cladding; h = {0}",h);
289
290 ileft = 0;
291 iright = points->axis[0]->size();
292 added_top_cladding = true;
293 }
294 }
295 }
296 }
297 if (!regions.empty() && regions.back().isQW(regions.back().size() - 1))
298 throw Exception("{0}: Quantum-well at the edge of the structure.", this->getId());
299
300 this->writelog(LOG_INFO, "Found {0} active region{1}", regions.size(), (regions.size() == 1) ? "" : "s");
301 size_t n = 0;
302 for (auto& region : regions) {
303 if (region.layers->getChildrenCount() <= 2) throw Exception("not enough layers in the active region {}", n);
304 region.summarize(this);
305 this->writelog(LOG_INFO, "Active region {0}: {1} nm single QW, {2} nm all QW, {3} nm total", n++, 0.1 * region.qwlen,
306 0.1 * region.qwtotallen, 0.1 * region.totallen);
307 }
308
309 // energy levels for active region with two identical QWs won't be calculated so QW widths must be changed
310 this->writelog(LOG_DETAIL, "Updating QW widths");
311 n = 0;
312 for (auto& region : regions) {
313 region.lens.clear();
314 int N = region.size(); // number of all layers in the active region (QW, barr, external)
315 int nQW = 0; // number of QWs counter
316 for (int i = 0; i < N; ++i) {
317 if (region.isQW(i)) {
318 nQW++;
319 region.QWs.insert(i - 1);
320 }
321 region.lens.push_back(region.getLayerBox(i).height() * 1e4); // in [A]
322 this->writelog(LOG_DEBUG, "Layer {0} thickness: {1} nm", i + 1, 0.1 * region.lens[i]);
323 }
324 this->writelog(LOG_DEBUG, "Number of QWs in the above active region: {0}", nQW);
325 this->writelog(LOG_DEBUG, "QW initial thickness: {0} nm", 0.1 * region.qwlen);
326
327 if (adjust_widths) {
328 double hstep = region.qwlen / QWwidthMod;
329 if (!(nQW % 2)) {
330 double dh0 = -(floor(nQW / 2) + 0.5) * hstep;
331 for (int i = 0; i < N; ++i) {
332 if (region.isQW(i)) {
333 region.lens[i] += dh0;
334 this->writelog(LOG_DEBUG, "Layer {0} thickness: {1} nm", i + 1, 0.1 * region.lens[i]);
335 dh0 += hstep;
336 }
337 }
338 } else {
339 double dh0 = -(int(nQW / 2)) * hstep;
340 for (int i = 0; i < N; ++i) {
341 if (region.isQW(i)) {
342 region.lens[i] += dh0;
343 this->writelog(LOG_DEBUG, "Layer {0} modified thickness: {1} nm", i + 1, 0.1 * region.lens[i]);
344 dh0 += hstep;
345 }
346 }
347 }
348 this->writelog(LOG_DEBUG, "QW thickness step: {0} nm", 0.1 * hstep);
349 }
350 n++;
351 }
352
353 return regions;
354}
355
356 template <typename GeometryType> void FermiNewGainSolver<GeometryType>::prepareActiveRegionsInfo() {
357 auto regs = detectActiveRegions(this->geometry);
358 regions.clear();
359 regions.reserve(regs.size());
360 regions.assign(regs.begin(), regs.end());
361 if (geometry_mod) {
362 regs = detectActiveRegions(this->geometry_mod);
363 if (regs.size() != regions.size())
364 throw Exception("modified geometry has different number of active regions ({}) than the main one ({})", regs.size(),
365 regions.size());
366 auto region = regions.begin();
367 for (const auto& reg : regs) (region++)->mod.reset(std::move(reg));
368 }
369}
370
371 template <typename GeometryType>
372 void FermiNewGainSolver<GeometryType>::findEnergyLevels(Levels& levels, const ActiveRegionInfo& region, double T) {
373 if (isnan(T) || T < 0.) throw ComputationError(this->getId(), "wrong temperature ({0}K)", T);
374
375 this->writelog(LOG_DETAIL, "Determining energy levels");
376
377 buildStructure(T, region, levels.bandsEc, levels.bandsEvhh, levels.bandsEvlh);
378
379 std::vector<double> dso;
380 dso.reserve(region.size());
381 for (int i = 0; i < region.size(); ++i) dso.push_back(region.getLayerMaterial(i)->Dso(T));
382
383 std::vector<kubly::struktura*> holes;
384 holes.reserve(2);
385 if (levels.bandsEvhh) holes.push_back(levels.bandsEvhh.get());
386 if (levels.bandsEvlh) holes.push_back(levels.bandsEvlh.get());
387
388 levels.Eg = region.getLayerMaterial(0)->CB(T, 0.) - region.getLayerMaterial(0)->VB(T, 0.); // cladding Eg (eV)
389
390 if (region.mod) {
391 buildStructure(T, *region.mod, levels.modbandsEc, levels.modbandsEvhh, levels.modbandsEvlh);
392
393 std::vector<kubly::struktura*> modholes;
394 modholes.reserve(2);
395 if (levels.modbandsEvhh) modholes.push_back(levels.modbandsEvhh.get());
396 if (levels.modbandsEvlh) modholes.push_back(levels.modbandsEvlh.get());
397 levels.activeRegion = plask::make_shared<kubly::obszar_aktywny>(levels.bandsEc.get(), holes, levels.modbandsEc.get(),
398 modholes, levels.Eg, dso, roughness, matrixElem, T);
399 } else {
400 levels.activeRegion =
401 plask::make_shared<kubly::obszar_aktywny>(levels.bandsEc.get(), holes, levels.Eg, dso, roughness, matrixElem, T);
402 }
403 levels.activeRegion->zrob_macierze_przejsc();
404}
405
406 template <typename GeometryType>
407 void FermiNewGainSolver<GeometryType>::buildStructure(double T,
408 const ActiveRegionData& region,
409 std::unique_ptr<kubly::struktura>& bandsEc,
410 std::unique_ptr<kubly::struktura>& bandsEvhh,
411 std::unique_ptr<kubly::struktura>& bandsEvlh) {
412 if (strains) {
413 if (!this->substrateMaterial) throw ComputationError(this->getId(), "no layer with role 'substrate' has been found");
414 if (maxLoglevel >= LOG_DEBUG) {
415 for (int i = 0; i < region.size(); ++i) {
416 double e = (this->substrateMaterial->lattC(T, 'a') - region.getLayerMaterial(i)->lattC(T, 'a')) /
417 region.getLayerMaterial(i)->lattC(T, 'a');
418 if ((i == 0) || (i == region.size() - 1)) e = 0.;
419 this->writelog(LOG_DEBUG, "Layer {0} - strain: {1}{2}", i + 1, e * 100., '%');
420 }
421 }
422 }
423
424 kubly::struktura* Ec = buildEc(T, region);
425 kubly::struktura* Evhh = buildEvhh(T, region);
426 kubly::struktura* Evlh = buildEvlh(T, region);
427
428 if (!Ec)
429 throw BadInput(this->getId(), "Conduction QW depth negative for electrons, check CB values of active-region materials");
430 if (!Evhh && !Evlh)
431 throw BadInput(this->getId(),
432 "Valence QW depth negative for both heavy holes and light holes, check VB values of "
433 "active-region materials");
434
435 bandsEc.reset(Ec);
436 bandsEc->gwiazdki.reserve(region.QWs.size());
437 bandsEc->gwiazdki.assign(region.QWs.begin(), region.QWs.end());
438
439 if (Evhh) {
440 bandsEvhh.reset(Evhh);
441 bandsEvhh->gwiazdki.reserve(region.QWs.size());
442 bandsEvhh->gwiazdki.assign(region.QWs.begin(), region.QWs.end());
443 }
444 if (Evlh) {
445 bandsEvlh.reset(Evlh);
446 bandsEvlh->gwiazdki.reserve(region.QWs.size());
447 bandsEvlh->gwiazdki.assign(region.QWs.begin(), region.QWs.end());
448 }
449}
450
451 template <typename GeometryType>
452 kubly::struktura* FermiNewGainSolver<GeometryType>::buildEc(double T, const ActiveRegionData& region) {
453 ptrVector<kubly::warstwa> levels;
454
455 int N = region.size(); // number of all layers in the active region (QW, barr, external)
456
457 double lattSub, straine = 0.;
458
459 if (strains) {
460 lattSub = this->substrateMaterial->lattC(T, 'a');
461 }
462
463 double DEc = region.getLayerMaterial(0)->CB(T, 0.); // Ec0 for cladding
464
465 double x = 0.;
466 double Ec = region.getLayerMaterial(0)->CB(T, 0.) - DEc;
467 this->writelog(LOG_DEBUG, "Layer {0} CB: {1} eV", 1, region.getLayerMaterial(0)->CB(T, 0.));
468 levels.data.emplace_back(new kubly::warstwa_skraj(kubly::warstwa_skraj::lewa, region.getLayerMaterial(0)->Me(T, 0.).c11,
469 region.getLayerMaterial(0)->Me(T, 0.).c00, x,
470 Ec)); // left cladding
471 for (int i = 1; i < N - 1; ++i) {
472 if (strains) straine = lattSub / region.getLayerMaterial(i)->lattC(T, 'a') - 1.;
473 double h = region.lens[i]; // tH (A)
474 double CBshift = 0.;
475 if (region.isQW(i)) CBshift = condQWshift;
476 Ec = region.getLayerMaterial(i)->CB(T, straine) + CBshift - DEc;
477 this->writelog(LOG_DEBUG, "Layer {0} CB: {1} eV", i + 1, region.getLayerMaterial(i)->CB(T, straine) + CBshift);
478 levels.data.emplace_back(new kubly::warstwa(region.getLayerMaterial(i)->Me(T, straine).c11,
479 region.getLayerMaterial(i)->Me(T, straine).c00, x, Ec, (x + h),
480 Ec)); // wells and barriers
481 x += h;
482 if (region.getLayerMaterial(i)->CB(T, straine) > DEc) return nullptr;
483 }
484
485 Ec = (region.getLayerMaterial(N - 1)->CB(T, 0.) - DEc);
486 this->writelog(LOG_DEBUG, "Layer {0} CB: {1} eV", N, region.getLayerMaterial(N - 1)->CB(T, 0.));
487 levels.data.emplace_back(new kubly::warstwa_skraj(kubly::warstwa_skraj::prawa, region.getLayerMaterial(N - 1)->Me(T, 0.).c11,
488 region.getLayerMaterial(N - 1)->Me(T, 0.).c00, x,
489 Ec)); // right cladding
490
491 this->writelog(LOG_DETAIL, "Computing energy levels for electrons");
492 return new kubly::struktura(levels.data, kubly::struktura::el);
493}
494
495 template <typename GeometryType>
496 kubly::struktura* FermiNewGainSolver<GeometryType>::buildEvhh(double T, const ActiveRegionData& region) {
497 ptrVector<kubly::warstwa> levels;
498
499 int N = region.size(); // number of all layers int the active region (QW, barr, external)
500
501 double lattSub, straine = 0.;
502
503 if (strains) {
504 lattSub = this->substrateMaterial->lattC(T, 'a');
505 }
506
507 double DEvhh = region.getLayerMaterial(0)->VB(T, 0., '*', 'H'); // Ev0 for cladding
508
509 double x = 0.;
510 double Evhh = -(region.getLayerMaterial(0)->VB(T, 0., '*', 'H') - DEvhh);
511 this->writelog(LOG_DEBUG, "Layer {0} VB(hh): {1} eV", 1, region.getLayerMaterial(0)->VB(T, 0., '*', 'H'));
512 levels.data.emplace_back(new kubly::warstwa_skraj(kubly::warstwa_skraj::lewa, region.getLayerMaterial(0)->Mhh(T, 0.).c11,
513 region.getLayerMaterial(0)->Mhh(T, 0.).c00, x,
514 Evhh)); // left cladding
515 for (int i = 1; i < N - 1; ++i) {
516 if (strains) straine = lattSub / region.getLayerMaterial(i)->lattC(T, 'a') - 1.;
517 double h = region.lens[i]; // h (A)
518 double VBshift = 0.;
519 if (region.isQW(i)) VBshift = valeQWshift;
520 Evhh = -(region.getLayerMaterial(i)->VB(T, straine, '*', 'H') + VBshift - DEvhh);
521 this->writelog(LOG_DEBUG, "Layer {0} VB(hh): {1} eV", i + 1,
522 region.getLayerMaterial(i)->VB(T, straine, '*', 'H') + VBshift);
523 levels.data.emplace_back(new kubly::warstwa(region.getLayerMaterial(i)->Mhh(T, straine).c11,
524 region.getLayerMaterial(i)->Mhh(T, straine).c00, x, Evhh, (x + h),
525 Evhh)); // wells and barriers
526 x += h;
527 if (region.getLayerMaterial(i)->VB(T, straine, '*', 'H') < DEvhh) return nullptr;
528 }
529
530 Evhh = -(region.getLayerMaterial(N - 1)->VB(T, 0., '*', 'H') - DEvhh);
531 this->writelog(LOG_DEBUG, "Layer {0} VB(hh): {1} eV", N, region.getLayerMaterial(N - 1)->VB(T, 0., '*', 'H'));
532 levels.data.emplace_back(new kubly::warstwa_skraj(kubly::warstwa_skraj::prawa, region.getLayerMaterial(N - 1)->Mhh(T, 0.).c11,
533 region.getLayerMaterial(N - 1)->Mhh(T, 0.).c00, x, Evhh));
534
535 this->writelog(LOG_DETAIL, "Computing energy levels for heavy holes");
536 return new kubly::struktura(levels.data, kubly::struktura::hh);
537}
538
539 template <typename GeometryType>
540 kubly::struktura* FermiNewGainSolver<GeometryType>::buildEvlh(double T, const ActiveRegionData& region) {
541 ptrVector<kubly::warstwa> levels;
542
543 int N = region.size(); // number of all layers int the active region (QW, barr, external)
544
545 double lattSub, straine = 0.;
546
547 if (strains) {
548 lattSub = this->substrateMaterial->lattC(T, 'a');
549 }
550
551 double DEvlh = region.getLayerMaterial(0)->VB(T, 0., '*', 'L'); // Ev0 for cladding
552
553 double x = 0.;
554 double Evlh = -(region.getLayerMaterial(0)->VB(T, 0., '*', 'L') - DEvlh);
555 this->writelog(LOG_DEBUG, "Layer {0} VB(lh): {1} eV", 1, region.getLayerMaterial(0)->VB(T, 0., '*', 'L'));
556 levels.data.emplace_back(new kubly::warstwa_skraj(kubly::warstwa_skraj::lewa, region.getLayerMaterial(0)->Mlh(T, 0.).c11,
557 region.getLayerMaterial(0)->Mlh(T, 0.).c00, x,
558 Evlh)); // left cladding
559 for (int i = 1; i < N - 1; ++i) {
560 if (strains) straine = lattSub / region.getLayerMaterial(i)->lattC(T, 'a') - 1.;
561 double h = region.lens[i]; // tH (A)
562 double VBshift = 0.;
563 if (region.isQW(i)) VBshift = valeQWshift;
564 Evlh = -(region.getLayerMaterial(i)->VB(T, straine, '*', 'L') + VBshift - DEvlh);
565 this->writelog(LOG_DEBUG, "Layer {0} VB(lh): {1} eV", i + 1,
566 region.getLayerMaterial(i)->VB(T, straine, '*', 'L') + VBshift);
567 levels.data.emplace_back(new kubly::warstwa(region.getLayerMaterial(i)->Mlh(T, straine).c11,
568 region.getLayerMaterial(i)->Mlh(T, straine).c00, x, Evlh, (x + h),
569 Evlh)); // wells and barriers
570 x += h;
571 if (region.getLayerMaterial(i)->VB(T, straine, '*', 'L') < DEvlh) return nullptr;
572 ;
573 }
574
575 Evlh = -(region.getLayerMaterial(N - 1)->VB(T, 0., '*', 'L') - DEvlh);
576 this->writelog(LOG_DEBUG, "Layer {0} VB(lh): {1} eV", N, region.getLayerMaterial(N - 1)->VB(T, 0., '*', 'L'));
577 levels.data.emplace_back(new kubly::warstwa_skraj(kubly::warstwa_skraj::prawa, region.getLayerMaterial(N - 1)->Mlh(T, 0.).c11,
578 region.getLayerMaterial(N - 1)->Mlh(T, 0.).c00, x, Evlh));
579
580 this->writelog(LOG_DETAIL, "Computing energy levels for light holes");
581 return new kubly::struktura(levels.data, kubly::struktura::lh);
582}
583
584 template <typename GeometryType>
585 void FermiNewGainSolver<GeometryType>::showEnergyLevels(std::string str,
586 const std::unique_ptr<kubly::struktura>& structure,
587 double nQW) {
588 std::vector<kubly::stan>::const_iterator it_stan = structure->rozwiazania.begin();
589 for (int iQW = 1; it_stan != structure->rozwiazania.end(); iQW++) {
590 bool calc_avg = true;
591 double avg_energy_level = 0.;
592 for (int i = 0; i < nQW; i++) {
593 avg_energy_level += it_stan->poziom;
594 this->writelog(plask::LOG_DETAIL, "QW {0} - energy level for {1}: {2} eV from cladding band edge", iQW, str,
595 it_stan->poziom);
596 it_stan++;
597 if (it_stan == structure->rozwiazania.end()) {
598 calc_avg = false;
599 break;
600 }
601 }
602 if (calc_avg)
603 this->writelog(plask::LOG_DETAIL, "QW {0} - average energy level for {1}: {2} eV from cladding band edge", iQW, str,
604 avg_energy_level / nQW);
605 }
606}
607
608 template <typename GeometryType>
609 kubly::wzmocnienie FermiNewGainSolver<GeometryType>::getGainModule(double wavelength,
610 double T,
611 double n,
612 const ActiveRegionInfo& region,
613 const Levels& levels) {
614 if (isnan(T) || T < 0.) throw ComputationError(this->getId(), "wrong temperature ({0}K)", T);
615 if (isnan(n)) throw ComputationError(this->getId(), "wrong carriers concentration ({0}/cm3)", n);
616 n = max(n, 1e-6); // To avoid hangs
617
618 double QWh = region.qwtotallen; // total thickness of QWs (Å)
619
620 // calculating nR
621 double QWnr = 0.;
622 int nQW = 0;
623 int N = region.size(); // number of all layers int the active region (QW, barr, external)
624 for (int i = 1; i < N - 1; ++i) {
625 if (region.isQW(i)) {
626 QWnr += region.getLayerMaterial(i)->nr(wavelength, T);
627 nQW++;
628 }
629 }
630 QWnr /= nQW;
631
632 double Eg = region.getLayerMaterial(0)->CB(T, 0.) - region.getLayerMaterial(0)->VB(T, 0.);
633 // TODO Add strain
634 double deltaEg = Eg - levels.Eg;
635
636 // this->writelog(LOG_DEBUG, "Creating gain module");
637 kubly::wzmocnienie gainModule(levels.activeRegion.get(), n * (QWh * 1e-8), T, QWnr, deltaEg, QWh, // QWh in Å
638 region.mod ? kubly::wzmocnienie::Z_POSZERZENIEM : kubly::wzmocnienie::Z_CHROPOWATOSCIA);
639
640 // this->writelog(LOG_DEBUG, "Recalculating carrier concentrations..");
641 // double step = 1e-1 * n; // TODO
642 // double concQW = n;
643 // double in1 = 0., tn1 = 0.;
644
645 // for (int i = 0; i < 5; i++) {
646 // while (1) {
647 // in1 = n * 10.;
648
649 // gainModule.nosniki_c = gainModule.nosniki_v = gainModule.przel_gest_z_cm2(in1 * QWh * 1e-8);
650 // double qFc = gainModule.qFlc;
651 // // double tFlv = gainModule.qFlv;
652
653 // std::vector<double> tN = levels.bandsEc->koncentracje_w_warstwach(qFc, T);
654 // auto maxElem = std::max_element(tN.begin(), tN.end());
655 // tn1 = *maxElem;
656 // tn1 = kubly::struktura::koncentracja_na_cm_3(tn1);
657 // if (tn1 >= concQW)
658 // n -= step;
659 // else {
660 // n += step;
661 // break;
662 // }
663 // }
664 // step *= 0.1;
665 // n -= step;
666 // }
667
668 if (maxLoglevel >= LOG_DEBUG) {
669 double qFc = gainModule.qFlc;
670 double qFv = gainModule.qFlv;
671 this->writelog(LOG_DEBUG, "Quasi-Fermi level for electrons: {0} eV from cladding conduction band edge", qFc);
672 this->writelog(LOG_DEBUG, "Quasi-Fermi level for holes: {0} eV from cladding valence band edge", -qFv);
673 std::vector<double> ne = levels.bandsEc->koncentracje_w_warstwach(qFc, T);
674 std::vector<double> nlh = levels.bandsEvlh->koncentracje_w_warstwach(-qFv, T);
675 std::vector<double> nhh = levels.bandsEvhh->koncentracje_w_warstwach(-qFv, T);
676 for (int i = 0; i <= (int)ne.size() - 1; i++)
677 this->writelog(LOG_DEBUG, "Carriers concentration in layer {:d} [cm(-3)]: el:{:.3e} lh:{:.3e} hh:{:.3e} ", i + 1,
678 kubly::struktura::koncentracja_na_cm_3(ne[i]), kubly::struktura::koncentracja_na_cm_3(nlh[i]),
679 kubly::struktura::koncentracja_na_cm_3(nhh[i]));
680 }
681
682 return gainModule;
683}
684
685 static const shared_ptr<OrderedAxis> zero_axis(new OrderedAxis({0.}));
686
688 template <typename GeometryT, typename T> struct DataBase : public LazyDataImpl<T> {
689 struct AveragedData {
690 shared_ptr<const RectangularMesh<2>> mesh;
691 LazyData<double> data;
692 double factor;
693 const FermiNewGainSolver<GeometryT>* solver;
694 const char* name;
695 AveragedData(const FermiNewGainSolver<GeometryT>* solver,
696 const char* name,
697 const shared_ptr<const MeshAxis>& haxis,
698 const typename FermiNewGainSolver<GeometryT>::ActiveRegionInfo& region)
699 : solver(solver), name(name) {
700 auto vaxis = plask::make_shared<OrderedAxis>();
701 for (size_t n = 0; n != region.size(); ++n) {
702 if (region.isQW(n)) {
703 auto box = region.getLayerBox(n);
704 vaxis->addPoint(0.5 * (box.lower.c1 + box.upper.c1));
705 }
706 }
707 mesh = plask::make_shared<const RectangularMesh<2>>(const_pointer_cast<MeshAxis>(haxis), vaxis,
708 RectangularMesh<2>::ORDER_01);
709 factor = 1. / vaxis->size();
710 }
711 size_t size() const { return mesh->axis[0]->size(); }
712 double operator[](size_t i) const {
713 double val = 0.;
714 for (size_t j = 0; j != mesh->axis[1]->size(); ++j) {
715 auto v = data[mesh->index(i, j)];
716 if (isnan(v)) throw ComputationError(solver->getId(), "wrong {0} ({1}) at {2}", name, v, mesh->at(i, j));
717 v = max(v, 1e-6); // To avoid hangs
718 val += v;
719 }
720 return val * factor;
721 }
722 };
723
724 FermiNewGainSolver<GeometryT>* solver;
725 std::vector<shared_ptr<MeshAxis>> regpoints;
726 std::vector<LazyData<T>> data;
727 shared_ptr<const MeshD<2>> dest_mesh;
728
729 void setupFromAxis(const shared_ptr<MeshAxis>& axis) {
730 regpoints.reserve(solver->regions.size());
731 InterpolationFlags flags(solver->geometry);
732 for (size_t r = 0; r != solver->regions.size(); ++r) {
733 std::set<double> pts;
734 auto box = solver->regions[r].getBoundingBox();
735 double y = 0.5 * (box.lower.c1 + box.upper.c1);
736 for (double x : *axis) {
737 auto p = flags.wrap(vec(x, y));
738 if (solver->regions[r].contains(p)) pts.insert(p.c0);
739 }
740 auto msh = plask::make_shared<OrderedAxis>();
741 msh->addOrderedPoints(pts.begin(), pts.end(), pts.size());
742 regpoints.emplace_back(std::move(msh));
743 }
744 }
745
746 DataBase(FermiNewGainSolver<GeometryT>* solver, const shared_ptr<const MeshD<2>>& dst_mesh)
747 : solver(solver), dest_mesh(dst_mesh) {
748 // Create horizontal points lists
749 if (solver->mesh) {
750 setupFromAxis(solver->mesh);
751 } else if (auto rect_mesh = dynamic_pointer_cast<const RectangularMesh<2>>(dst_mesh)) {
752 setupFromAxis(rect_mesh->axis[0]);
753 } else {
754 regpoints.reserve(solver->regions.size());
755 InterpolationFlags flags(solver->geometry);
756 for (size_t r = 0; r != solver->regions.size(); ++r) {
757 std::set<double> pts;
758 for (auto point : *dest_mesh) {
759 auto p = flags.wrap(point);
760 if (solver->regions[r].contains(p)) pts.insert(p.c0);
761 }
762 auto msh = plask::make_shared<OrderedAxis>();
763 msh->addOrderedPoints(pts.begin(), pts.end(), pts.size());
764 regpoints.emplace_back(std::move(msh));
765 }
766 }
767 }
768
769 void compute(double wavelength, InterpolationMethod interp) {
770 // Compute gains on mesh for each active region
771 data.resize(solver->regions.size());
772 for (size_t reg = 0; reg != solver->regions.size(); ++reg) {
773 if (regpoints[reg]->size() == 0) {
774 data[reg] = LazyData<T>(dest_mesh->size(), Zero<T>());
775 continue;
776 }
777 DataVector<T> values(regpoints[reg]->size());
778 AveragedData temps(solver, "temperature", regpoints[reg], solver->regions[reg]);
779 AveragedData concs(temps);
780 concs.name = "carriers concentration";
781 temps.data = solver->inTemperature(temps.mesh, interp);
782 concs.data = solver->inCarriersConcentration(temps.mesh, interp);
783 if (solver->build_struct_once && !solver->region_levels[reg]) {
784 if (isnan(solver->Tref))
785 throw ComputationError(solver->getId(), "no reference temperature set for fast levels calculation");
786 solver->findEnergyLevels(solver->region_levels[reg], solver->regions[reg], solver->Tref);
787 }
788 std::exception_ptr error;
789 PLASK_OMP_PARALLEL_FOR
790 for (int i = 0; i < regpoints[reg]->size(); ++i) {
791 if (error) continue;
792 try {
793 if (solver->build_struct_once) {
794 values[i] =
795 getValue(wavelength, temps[i], max(concs[i], 1e-9), solver->regions[reg], solver->region_levels[reg]);
796 } else {
797 Levels levels;
798 solver->findEnergyLevels(levels, solver->regions[reg], temps[i]);
799 values[i] = getValue(wavelength, temps[i], max(concs[i], 1e-9), solver->regions[reg], levels);
800 }
801 } catch (...) {
802#pragma omp critical
803 error = std::current_exception();
804 }
805 }
806 if (error) std::rethrow_exception(error);
807 data[reg] = interpolate(plask::make_shared<RectangularMesh<2>>(regpoints[reg], zero_axis), values, dest_mesh, interp);
808 }
809 }
810
811 virtual T getValue(double wavelength,
812 double temp,
813 double conc,
814 const typename FermiNewGainSolver<GeometryT>::ActiveRegionInfo& region,
815 const Levels& levels) = 0;
816
817 size_t size() const override { return dest_mesh->size(); }
818
819 T at(size_t i) const override {
820 for (size_t reg = 0; reg != solver->regions.size(); ++reg)
821 // if (solver->regions[reg].contains(dest_mesh->at(i)))
822 if (solver->regions[reg].inQW(dest_mesh->at(i))) return data[reg][i];
823 return Zero<T>();
824 }
825};
826
827 template <typename GeometryT> struct GainData : public DataBase<GeometryT, Tensor2<double>> {
828 template <typename... Args> GainData(Args... args) : DataBase<GeometryT, Tensor2<double>>(args...) {}
829
830 Tensor2<double> getValue(double wavelength,
831 double temp,
832 double conc,
833 const typename FermiNewGainSolver<GeometryT>::ActiveRegionInfo& region,
834 const Levels& levels) override {
835 kubly::wzmocnienie gainModule = this->solver->getGainModule(wavelength, temp, conc, region, levels);
836
837 double E = nm_to_eV(wavelength);
838 if (!this->solver->lifetime || region.mod)
839 return Tensor2<double>(gainModule.wzmocnienie_calk_bez_splotu(E, 0.), gainModule.wzmocnienie_calk_bez_splotu(E, 1.));
840 else {
841 double beta = phys::hb_eV * 1e12 / this->solver->lifetime;
842 return Tensor2<double>(gainModule.wzmocnienie_calk_ze_splotem(E, beta, 0.),
843 gainModule.wzmocnienie_calk_ze_splotem(E, beta, 1.));
844 }
845 }
846};
847
848 template <typename GeometryT> struct DgDnData : public DataBase<GeometryT, Tensor2<double>> {
849 template <typename... Args> DgDnData(Args... args) : DataBase<GeometryT, Tensor2<double>>(args...) {}
850
851 Tensor2<double> getValue(double wavelength,
852 double temp,
853 double conc,
854 const typename FermiNewGainSolver<GeometryT>::ActiveRegionInfo& region,
855 const Levels& levels) override {
856 double h = 0.5 * this->solver->differenceQuotient;
857 double conc1 = (1. - h) * conc, conc2 = (1. + h) * conc;
858
859 kubly::wzmocnienie gainModule1 = this->solver->getGainModule(wavelength, temp, conc1, region, levels);
860 kubly::wzmocnienie gainModule2 = this->solver->getGainModule(wavelength, temp, conc2, region, levels);
861
862 double E = nm_to_eV(wavelength);
863 Tensor2<double> gain1, gain2;
864 if (!this->solver->lifetime || region.mod) {
865 gain1 = Tensor2<double>(gainModule1.wzmocnienie_calk_bez_splotu(E, 0.), //
866 gainModule1.wzmocnienie_calk_bez_splotu(E, 1.));
867 gain2 = Tensor2<double>(gainModule2.wzmocnienie_calk_bez_splotu(E, 0.), //
868 gainModule2.wzmocnienie_calk_bez_splotu(E, 1.));
869 } else {
870 double beta = phys::hb_eV * 1e12 / this->solver->lifetime;
871 gain1 = Tensor2<double>(gainModule1.wzmocnienie_calk_ze_splotem(E, beta, 0.),
872 gainModule1.wzmocnienie_calk_ze_splotem(E, beta, 1.));
873 gain2 = Tensor2<double>(gainModule2.wzmocnienie_calk_ze_splotem(E, beta, 0.),
874 gainModule2.wzmocnienie_calk_ze_splotem(E, beta, 1.));
875 }
876 return (gain2 - gain1) / (2. * h * conc);
877 }
878};
879
880 inline static Tensor2<double> sumLuminescence(kubly::wzmocnienie& gain, double wavelength) {
881 double E = nm_to_eV(wavelength);
882 double resultTE = 0., resultTM = 0.;
883 for (int nr_c = 0; nr_c <= (int)gain.pasma->pasmo_przew.size() - 1; nr_c++) {
884 for (int nr_v = 0; nr_v <= (int)gain.pasma->pasmo_wal.size() - 1; nr_v++) {
885 resultTE += gain.spont_od_pary_pasm(E, nr_c, nr_v, 0.0);
886 resultTM += gain.spont_od_pary_pasm(E, nr_c, nr_v, 1.0);
887 }
888 }
889 return Tensor2<double>(resultTE, resultTM);
890}
891
892 template <typename GeometryT> struct LuminescenceData : public DataBase<GeometryT, Tensor2<double>> {
893 template <typename... Args> LuminescenceData(Args... args) : DataBase<GeometryT, Tensor2<double>>(args...) {}
894
895 Tensor2<double> getValue(double wavelength,
896 double temp,
897 double conc,
898 const typename FermiNewGainSolver<GeometryT>::ActiveRegionInfo& region,
899 const Levels& levels) override {
900 kubly::wzmocnienie gainModule = this->solver->getGainModule(wavelength, temp, conc, region, levels);
901
902 double QWfrac = region.qwtotallen / region.totallen;
903 return sumLuminescence(gainModule, wavelength) / QWfrac;
904 }
905};
906
907 template <typename GeometryType>
908 const LazyData<Tensor2<double>> FermiNewGainSolver<GeometryType>::getGain(Gain::EnumType what,
909 const shared_ptr<const MeshD<2>>& dst_mesh,
910 double wavelength,
911 InterpolationMethod interp) {
912 if (what == Gain::DGDN) {
913 this->writelog(LOG_DETAIL, "Calculating gain over carriers concentration derivative");
914 this->initCalculation(); // This must be called before any calculation!
915 DgDnData<GeometryType>* data = new DgDnData<GeometryType>(this, dst_mesh);
916 data->compute(wavelength, getInterpolationMethod<INTERPOLATION_SPLINE>(interp));
917 return LazyData<Tensor2<double>>(data);
918 } else {
919 this->writelog(LOG_DETAIL, "Calculating gain");
920 this->initCalculation(); // This must be called before any calculation!
921 GainData<GeometryType>* data = new GainData<GeometryType>(this, dst_mesh);
922 data->compute(wavelength, getInterpolationMethod<INTERPOLATION_SPLINE>(interp));
923 return LazyData<Tensor2<double>>(data);
924 }
925}
926
927 template <typename GeometryType>
928 const LazyData<Tensor2<double>> FermiNewGainSolver<GeometryType>::getLuminescence(const shared_ptr<const MeshD<2>>& dst_mesh,
929 double wavelength,
930 InterpolationMethod interp) {
931 this->writelog(LOG_DETAIL, "Calculating luminescence");
932 this->initCalculation(); // This must be called before any calculation!
933
934 LuminescenceData<GeometryType>* data = new LuminescenceData<GeometryType>(this, dst_mesh);
935 data->compute(wavelength, getInterpolationMethod<INTERPOLATION_SPLINE>(interp));
936
937 return LazyData<Tensor2<double>>(data);
938}
939
940 template <typename GeometryT>
941 GainSpectrum<GeometryT>::GainSpectrum(FermiNewGainSolver<GeometryT>* solver, const Vec<2> point) : solver(solver), point(point) {
942 auto mesh = plask::make_shared<const OnePointMesh<2>>(point);
943 T = solver->inTemperature(mesh)[0];
944 n = solver->inCarriersConcentration(mesh)[0];
945 for (reg = 0; reg < solver->regions.size(); ++reg) {
946 if (solver->regions[reg].contains(point)) {
947 solver->inTemperature.changedConnectMethod(this, &GainSpectrum::onTChange);
948 solver->inCarriersConcentration.changedConnectMethod(this, &GainSpectrum::onNChange);
949 return;
950 };
951 }
952 throw BadInput(solver->getId(), "point {0} does not belong to any active region", point);
953}
954
955 template <typename GeometryT> Tensor2<double> GainSpectrum<GeometryT>::getGain(double wavelength) {
956 if (!gMod) {
957 Levels* levels;
958 if (solver->build_struct_once) {
959 if (!solver->region_levels[reg]) {
960 if (isnan(solver->Tref))
961 throw ComputationError(solver->getId(), "no reference temperature set for fast levels calculation");
962 solver->findEnergyLevels(solver->region_levels[reg], solver->regions[reg], solver->Tref);
963 }
964 levels = &solver->region_levels[reg];
965 } else {
966 this->levels.reset(new Levels);
967 levels = this->levels.get();
968 solver->findEnergyLevels(*levels, solver->regions[reg], T);
969 }
970 gMod.reset(new kubly::wzmocnienie(std::move(solver->getGainModule(wavelength, T, n, solver->regions[reg], *levels))));
971 }
972
973 double E = nm_to_eV(wavelength);
974 double tau = solver->getLifeTime();
975 if (!tau || solver->regions[reg].mod)
976 return Tensor2<double>(gMod->wzmocnienie_calk_bez_splotu(E, 0.), gMod->wzmocnienie_calk_bez_splotu(E, 1.));
977 else {
978 double beta = phys::hb_eV * 1e12 / tau;
979 return Tensor2<double>(gMod->wzmocnienie_calk_ze_splotem(E, beta, 0.), gMod->wzmocnienie_calk_ze_splotem(E, beta, 1.));
980 }
981}
982
983 template <typename GeometryT>
984 LuminescenceSpectrum<GeometryT>::LuminescenceSpectrum(FermiNewGainSolver<GeometryT>* solver, const Vec<2> point)
985 : solver(solver), point(point) {
986 auto mesh = plask::make_shared<const OnePointMesh<2>>(point);
987 T = solver->inTemperature(mesh)[0];
988 n = solver->inCarriersConcentration(mesh)[0];
989 for (reg = 0; reg < solver->regions.size(); ++reg) {
990 if (solver->regions[reg].contains(point)) {
991 solver->inTemperature.changedConnectMethod(this, &LuminescenceSpectrum::onTChange);
992 solver->inCarriersConcentration.changedConnectMethod(this, &LuminescenceSpectrum::onNChange);
993 return;
994 };
995 }
996 throw BadInput(solver->getId(), "point {0} does not belong to any active region", point);
997}
998
999 template <typename GeometryT> Tensor2<double> LuminescenceSpectrum<GeometryT>::getLuminescence(double wavelength) {
1000 if (!gMod) {
1001 Levels* levels;
1002 if (solver->build_struct_once) {
1003 if (!solver->region_levels[reg]) {
1004 if (isnan(solver->Tref))
1005 throw ComputationError(solver->getId(), "no reference temperature set for fast levels calculation");
1006 solver->findEnergyLevels(solver->region_levels[reg], solver->regions[reg], solver->Tref);
1007 }
1008 levels = &solver->region_levels[reg];
1009 } else {
1010 this->levels.reset(new Levels);
1011 levels = this->levels.get();
1012 solver->findEnergyLevels(*levels, solver->regions[reg], T);
1013 }
1014 gMod.reset(new kubly::wzmocnienie(std::move(solver->getGainModule(wavelength, T, n, solver->regions[reg], *levels))));
1015 }
1016 double QWfrac = solver->regions[reg].qwtotallen / solver->regions[reg].totallen;
1017 return sumLuminescence(*gMod, wavelength) / QWfrac;
1018}
1019
1020 template <typename GeometryType> GainSpectrum<GeometryType> FermiNewGainSolver<GeometryType>::getGainSpectrum(const Vec<2>& point) {
1021 this->initCalculation();
1022 return GainSpectrum<GeometryType>(this, point);
1023}
1024
1025 template <typename GeometryType>
1026 LuminescenceSpectrum<GeometryType> FermiNewGainSolver<GeometryType>::getLuminescenceSpectrum(const Vec<2>& point) {
1027 this->initCalculation();
1028 return LuminescenceSpectrum<GeometryType>(this, point);
1029}
1030
1031 template <> std::string FermiNewGainSolver<Geometry2DCartesian>::getClassName() const { return "gain.FermiNew2D"; }
1032 template <> std::string FermiNewGainSolver<Geometry2DCylindrical>::getClassName() const { return "gain.FermiNewCyl"; }
1033
1034 template struct PLASK_SOLVER_API FermiNewGainSolver<Geometry2DCartesian>;
1035 template struct PLASK_SOLVER_API FermiNewGainSolver<Geometry2DCylindrical>;
1036
1037 template struct PLASK_SOLVER_API GainSpectrum<Geometry2DCartesian>;
1038 template struct PLASK_SOLVER_API GainSpectrum<Geometry2DCylindrical>;
1039
1040 template struct PLASK_SOLVER_API LuminescenceSpectrum<Geometry2DCartesian>;
1041 template struct PLASK_SOLVER_API LuminescenceSpectrum<Geometry2DCylindrical>;
1042
1043}}} // namespace plask::solvers::FermiNew