PLaSK library

Loading...
Searching...
No Matches
efm.hpp
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#ifndef PLASK__MODULE_OPTICAL_EFM_HPP
15#define PLASK__MODULE_OPTICAL_EFM_HPP
16
17#include <limits>
18
19#include <camos/camos.h>
20#include <plask/plask.hpp>
21
22#include "bisection.hpp"
23#include "rootdigger.hpp"
24
25 namespace plask { namespace optical { namespace effective {
26
27 static constexpr int MH = 2; // Hankel function type (1 or 2)
28
32 struct PLASK_SOLVER_API EffectiveFrequencyCyl : public SolverWithMesh<Geometry2DCylindrical, RectangularMesh<2>> {
33 struct FieldZ {
34 dcomplex F, B;
35 FieldZ() = default;
36 FieldZ(dcomplex f, dcomplex b) : F(f), B(b) {}
37 FieldZ operator*(dcomplex a) const { return FieldZ(F * a, B * a); }
38 FieldZ operator/(dcomplex a) const { return FieldZ(F / a, B / a); }
39 FieldZ operator*=(dcomplex a) {
40 F *= a;
41 B *= a;
42 return *this;
43 }
44 FieldZ operator/=(dcomplex a) {
45 F /= a;
46 B /= a;
47 return *this;
48 }
49 };
50
51 struct MatrixZ {
52 dcomplex ff, fb, bf, bb;
53 MatrixZ() = default;
54 MatrixZ(dcomplex t1, dcomplex t2, dcomplex t3, dcomplex t4) : ff(t1), fb(t2), bf(t3), bb(t4) {}
55 static MatrixZ eye() { return MatrixZ(1., 0., 0., 1.); }
56 static MatrixZ diag(dcomplex f, dcomplex b) { return MatrixZ(f, 0., 0., b); }
57 MatrixZ operator*(const MatrixZ& T) {
58 return MatrixZ(ff * T.ff + fb * T.bf, ff * T.fb + fb * T.bb, bf * T.ff + bb * T.bf, bf * T.fb + bb * T.bb);
59 }
60 FieldZ operator*(const FieldZ& v) { return FieldZ(ff * v.F + fb * v.B, bf * v.F + bb * v.B); }
61 FieldZ solve(const FieldZ& v) { return FieldZ(bb * v.F - fb * v.B, -bf * v.F + ff * v.B) / (ff * bb - fb * bf); }
62 };
63
64 struct FieldR {
65 dcomplex J, H;
66 FieldR() = default;
67 FieldR(dcomplex j, dcomplex h) : J(j), H(h) {}
68 FieldR operator*(dcomplex a) const { return FieldR(a * J, a * H); }
69 FieldR operator/(dcomplex a) const { return FieldR(J / a, H / a); }
70 FieldR& operator*=(dcomplex a) {
71 J *= a;
72 H *= a;
73 return *this;
74 }
75 FieldR& operator/=(dcomplex a) {
76 J /= a;
77 H /= a;
78 return *this;
79 }
80 };
81
82 struct MatrixR {
83 dcomplex JJ, JH, HJ, HH;
84 MatrixR(dcomplex jj, dcomplex jh, dcomplex hj, dcomplex hh) : JJ(jj), JH(jh), HJ(hj), HH(hh) {}
85 static MatrixR eye() { return MatrixR(1., 0., 0., 1.); }
86 static MatrixR diag(dcomplex j, dcomplex h) { return MatrixR(j, 0., 0., h); }
87 MatrixR operator*(dcomplex c) { return MatrixR(c * JJ, c * JH, c * HJ, c * HH); }
88 friend MatrixR operator*(dcomplex c, const MatrixR& M) { return MatrixR(c * M.JJ, c * M.JH, c * M.HJ, c * M.HH); }
89 MatrixR operator/(dcomplex d) {
90 dcomplex c = 1. / d;
91 return MatrixR(c * JJ, c * JH, c * HJ, c * HH);
92 }
93 MatrixR& operator*=(dcomplex c) {
94 JJ *= c;
95 JH *= c;
96 HJ *= c;
97 HH *= c;
98 return *this;
99 }
100 MatrixR& operator/=(dcomplex d) {
101 dcomplex c = 1. / d;
102 JJ *= c;
103 JH *= c;
104 HJ *= c;
105 HH *= c;
106 return *this;
107 }
108 FieldR operator*(const FieldR& v) { return FieldR(JJ * v.J + JH * v.H, HJ * v.J + HH * v.H); }
109 MatrixR operator*(const MatrixR& o) {
110 return MatrixR(JJ * o.JJ + JH * o.HJ, JJ * o.JH + JH * o.HH, HJ * o.JJ + HH * o.HJ, HJ * o.JH + HH * o.HH);
111 }
112 FieldR solve(const FieldR& v) { return FieldR(HH * v.J - JH * v.H, -HJ * v.J + JJ * v.H) / (JJ * HH - JH * HJ); }
113 MatrixR solve(const MatrixR& o) {
114 MatrixR result(HH * o.JJ - JH * o.HJ, HH * o.JH - JH * o.HH, -HJ * o.JJ + JJ * o.HJ, -HJ * o.JH + JJ * o.HH);
115 result /= (JJ * HH - JH * HJ);
116 return result;
117 }
118 };
119
121 enum Emission {
122 TOP,
123 BOTTOM
124 };
125
127 enum Determinant {
128 DETERMINANT_INWARDS,
129 DETERMINANT_OUTWARDS,
130 DETERMINANT_FULL
131 };
132
134 struct Mode {
135 EffectiveFrequencyCyl* solver;
136 int m;
137 bool have_fields;
138 std::vector<FieldR, aligned_allocator<FieldR>> rfields;
139 std::vector<double, aligned_allocator<double>> rweights;
140 dcomplex lam;
141 double power;
142
143 Mode(EffectiveFrequencyCyl* solver, int m = 0)
144 : solver(solver), m(m), have_fields(false), rfields(solver->rsize), rweights(solver->rsize), power(1.) {}
145
146 bool operator==(const Mode& other) const { return m == other.m && is_zero(lam - other.lam); }
147
149 dcomplex rField(double r) const {
150 double Jr, Ji, Hr, Hi;
151 long nz, ierr;
152 size_t ir = solver->mesh->axis[0]->findIndex(r);
153 if (ir > 0) --ir;
154 if (ir >= solver->veffs.size()) ir = solver->veffs.size() - 1;
155 dcomplex x = r * solver->k0 * sqrt(solver->nng[ir] * (solver->veffs[ir] - solver->freqv(lam)));
156 if (real(x) < 0.) x = -x;
157 if (imag(x) > SMALL) x = -x;
158 if (ir == solver->rsize - 1) {
159 Jr = Ji = 0.;
160 } else {
161 zbesj(x.real(), x.imag(), m, 1, 1, &Jr, &Ji, nz, ierr);
162 if (ierr != 0) throw ComputationError(solver->getId(), "could not compute J({0}, {1}) @ r = {2}um", m, str(x), r);
163 }
164 if (ir == 0) {
165 Hr = Hi = 0.;
166 } else {
167 zbesh(x.real(), x.imag(), m, 1, MH, 1, &Hr, &Hi, nz, ierr);
168 if (ierr != 0) throw ComputationError(solver->getId(), "could not compute H({0}, {1}) @ r = {2}um", m, str(x), r);
169 }
170 return rfields[ir].J * dcomplex(Jr, Ji) + rfields[ir].H * dcomplex(Hr, Hi);
171 }
172
174 double loss() const { return imag(4e7 * PI / lam); }
175 };
176
178 dcomplex freqv(dcomplex lam) { return 2. - 4e3 * PI / lam / k0; }
179
181 dcomplex lambda(dcomplex freq) { return 2e3 * PI / (k0 * (1. - freq / 2.)); }
182
183 protected:
184 friend struct RootDigger;
185
187 DataLog<dcomplex, dcomplex> log_value;
188
189 size_t rsize,
190 zbegin,
191 zsize;
192
194 std::vector<std::vector<dcomplex, aligned_allocator<dcomplex>>> nrCache;
195
197 std::vector<std::vector<dcomplex, aligned_allocator<dcomplex>>> ngCache;
198
200 std::vector<FieldZ> zfields;
201
203 std::vector<double, aligned_allocator<double>> zintegrals;
204
206 std::vector<dcomplex, aligned_allocator<dcomplex>> veffs;
207
209 std::vector<dcomplex, aligned_allocator<dcomplex>> nng;
210
212 dcomplex old_k0;
213
215 Emission emission;
216
218 void onInputChange(ReceiverBase&, ReceiverBase::ChangeReason) { cache_outdated = true; }
219
224 int rstripe;
225
226 public:
228 Determinant determinant;
229
231 int getStripe() const { return rstripe; }
232
234 void setStripe(int stripe) {
235 if (!mesh) setSimpleMesh();
236 if (stripe < 0 || std::size_t(stripe) >= mesh->axis[0]->size()) throw BadInput(getId(), "wrong stripe number specified");
237 rstripe = stripe;
238 invalidate();
239 }
240
242 double getStripeR() const {
243 if (rstripe == -1 || !mesh) return NAN;
244 return mesh->axis[0]->at(rstripe);
245 }
246
251 void setStripeR(double r = 0.) {
252 if (!mesh) setSimpleMesh();
253 if (r < 0) throw BadInput(getId(), "radial position cannot be negative");
254 rstripe = int(std::lower_bound(mesh->axis[0]->begin() + 1, mesh->axis[0]->end(), r) - mesh->axis[0]->begin() - 1);
255 invalidate();
256 }
257
259 void useAllStripes() {
260 rstripe = -1;
261 invalidate();
262 }
263
268 dcomplex getDeltaNeff(double r) {
269 stageOne();
270 if (r < 0) throw BadInput(getId(), "radial position cannot be negative");
271 size_t ir = mesh->axis[0]->findIndex(r);
272 if (ir > 0) --ir;
273 if (ir >= veffs.size()) ir = veffs.size() - 1;
274 return sqrt(nng[ir] * veffs[ir]);
275 }
276
281 dcomplex getNNg(double r) {
282 stageOne();
283 if (r < 0) throw BadInput(getId(), "radial position cannot be negative");
284 size_t ir = mesh->axis[0]->findIndex(r);
285 if (ir > 0) --ir;
286 if (ir >= veffs.size()) ir = veffs.size() - 1;
287 return sqrt(nng[ir]);
288 }
289
290 // Parameters for rootdigger
291 RootDigger::Params root;
292 RootDigger::Params stripe_root;
293
295 double perr;
296
298 dcomplex k0;
299
301 dcomplex vlam;
302
304 std::vector<Mode> modes;
305
307 ReceiverFor<Temperature, Geometry2DCylindrical> inTemperature;
308
310 ReceiverFor<Gain, Geometry2DCylindrical> inGain;
311
313 ReceiverFor<CarriersConcentration, Geometry2DCylindrical> inCarriersConcentration;
314
316 typename ProviderFor<ModeWavelength>::Delegate outWavelength;
317
319 typename ProviderFor<ModeLoss>::Delegate outLoss;
320
322 typename ProviderFor<ModeLightMagnitude, Geometry2DCylindrical>::Delegate outLightMagnitude;
323
325 typename ProviderFor<ModeLightE, Geometry2DCylindrical>::Delegate outLightE;
326
328 typename ProviderFor<RefractiveIndex, Geometry2DCylindrical>::Delegate outRefractiveIndex;
329
331 typename ProviderFor<Heat, Geometry2DCylindrical>::Delegate outHeat;
332
333 EffectiveFrequencyCyl(const std::string& name = "");
334
335 virtual ~EffectiveFrequencyCyl() {
336 inTemperature.changedDisconnectMethod(this, &EffectiveFrequencyCyl::onInputChange);
337 inGain.changedDisconnectMethod(this, &EffectiveFrequencyCyl::onInputChange);
338 inCarriersConcentration.changedDisconnectMethod(this, &EffectiveFrequencyCyl::onInputChange);
339 }
340
341 std::string getClassName() const override { return "optical.EffectiveFrequencyCyl"; }
342
343 std::string getClassDescription() const override {
344 return "Calculate optical modes and optical field distribution using the effective index method "
345 "in Cartesian two-dimensional space.";
346 }
347
348 void loadConfiguration(plask::XMLReader& reader, plask::Manager& manager) override;
349
352 Emission getEmission() const { return emission; }
353
356 void setEmission(Emission emis) {
357 emission = emis;
358 for (auto& mode : modes) mode.have_fields = false;
359 }
360
364 void setSimpleMesh() {
365 writelog(LOG_DETAIL, "Creating simple mesh");
366 setMesh(plask::make_shared<RectangularMesh2DSimpleGenerator>());
367 }
368
375 void setHorizontalMesh(shared_ptr<MeshAxis> meshx) {
376 writelog(LOG_DETAIL, "Setting horizontal mesh");
377 if (!geometry) throw NoChildException();
378 auto meshxy = makeGeometryGrid(geometry->getChild());
379 meshxy->setTran(meshx);
380 setMesh(meshxy);
381 }
382
384 bool getAsymptotic() const { return asymptotic; }
385
387 void setAsymptotic(bool value) {
388 asymptotic = value;
389 invalidate();
390 }
391
401 size_t findMode(dcomplex lambda, int m = 0);
402
416 std::vector<size_t> findModes(plask::dcomplex lambda1 = 0.,
417 plask::dcomplex lambda2 = 0.,
418 int m = 0,
419 size_t resteps = 256,
420 size_t imsteps = 64,
421 dcomplex eps = dcomplex(1e-6, 1e-9));
422
427 dcomplex getVertDeterminant(dcomplex vlambda) {
428 updateCache();
429 if (rstripe < 0) throw BadInput(getId(), "this works only for the weighted approach");
430 if (vlam == 0. && isnan(k0.real())) throw BadInput(getId(), "no reference wavelength `lam0` specified");
431 dcomplex v = freqv(vlambda);
432 return this->detS1(v, nrCache[rstripe], ngCache[rstripe]);
433 }
434
440 dcomplex getDeterminant(dcomplex lambda, int m = 0) {
441 if (isnan(k0.real())) throw BadInput(getId(), "no reference wavelength `lam0` specified");
442 stageOne();
443 Mode mode(this, m);
444 dcomplex det = detS(lambda, mode);
445 // log_value(v, det);
446 return det;
447 }
448
455 size_t setMode(dcomplex clambda, int m = 0);
456
465 inline size_t setMode(double lambda, double loss, int m = 0) {
466 return setMode(dcomplex(lambda, -lambda * lambda / (4e7 * PI) * loss), m);
467 }
468
470 void clearModes() { modes.clear(); }
471
476 double getTotalAbsorption(Mode& mode);
477
482 double getTotalAbsorption(size_t num);
483
488 double getGainIntegral(Mode& mode);
489
494 double getGainIntegral(size_t num);
495
496 protected:
498 bool need_gain;
499
501 bool cache_outdated;
502
504 bool have_veffs;
505
507 bool asymptotic;
508
510 void onInitialize() override;
511
513 void onInvalidate() override;
514
518 void updateCache();
519
524 void stageOne();
525
527 dcomplex detS1(const dcomplex& v,
528 const std::vector<dcomplex, aligned_allocator<dcomplex>>& NR,
529 const std::vector<dcomplex, aligned_allocator<dcomplex>>& NG,
530 std::vector<FieldZ>* saveto = nullptr);
531
533 void computeStripeNNg(size_t stripe, bool save_integrals = false);
534
538 double integrateBessel(Mode& mode);
539
541 void computeBessel(size_t i, dcomplex v, const Mode& mode, dcomplex* J1, dcomplex* H1, dcomplex* J2, dcomplex* H2);
542
544 dcomplex detS(const plask::dcomplex& lam, Mode& mode, bool save = false);
545
547 size_t getMainStripe() {
548 if (rstripe < 0) {
549 size_t stripe = 0;
550 // Look for the innermost stripe with not constant refractive index
551 bool all_the_same = true;
552 while (all_the_same) {
553 dcomplex same_nr = nrCache[stripe].front();
554 dcomplex same_ng = ngCache[stripe].front();
555 for (auto nr = nrCache[stripe].begin(), ng = ngCache[stripe].begin(); nr != nrCache[stripe].end(); ++nr, ++ng)
556 if (*nr != same_nr || *ng != same_ng) {
557 all_the_same = false;
558 break;
559 }
560 if (all_the_same) ++stripe;
561 }
562 writelog(LOG_DETAIL, "Vertical field distribution taken from stripe {0}", stripe);
563 return stripe;
564 } else {
565 return rstripe;
566 }
567 }
568
570 size_t insertMode(const Mode& mode) {
571 for (size_t i = 0; i != modes.size(); ++i)
572 if (modes[i] == mode) return i;
573 modes.push_back(mode);
574 outWavelength.fireChanged();
575 outLoss.fireChanged();
576 outLightMagnitude.fireChanged();
577 outLightE.fireChanged();
578 return modes.size() - 1;
579 }
580
582 size_t nmodes() const { return modes.size(); }
583
588 double getWavelength(size_t n) {
589 if (n >= modes.size()) throw NoValue(ModeWavelength::NAME);
590 return real(modes[n].lam);
591 }
592
597 double getModalLoss(size_t n) {
598 if (n >= modes.size()) throw NoValue(ModeLoss::NAME);
599 return imag(4e7 * PI / modes[n].lam); // 2e4 2/µm -> 2/cm
600 }
601
602 template <typename T> struct FieldDataBase;
603 template <typename T> struct FieldDataInefficient;
604 template <typename T> struct FieldDataEfficient;
605 struct HeatDataImpl;
606
608 const LazyData<double> getLightMagnitude(std::size_t num,
609 const shared_ptr<const MeshD<2>>& dst_mesh,
610 InterpolationMethod = INTERPOLATION_DEFAULT);
611
613 const LazyData<Vec<3, dcomplex>> getElectricField(std::size_t num,
614 const shared_ptr<const plask::MeshD<2>>& dst_mesh,
615 plask::InterpolationMethod = INTERPOLATION_DEFAULT);
616
618 const LazyData<dcomplex> getRefractiveIndex(RefractiveIndex::EnumType component,
619 const shared_ptr<const MeshD<2>>& dst_mesh,
620 dcomplex lam,
621 InterpolationMethod = INTERPOLATION_DEFAULT);
622
624 const LazyData<double> getHeat(const shared_ptr<const MeshD<2>>& dst_mesh, InterpolationMethod method = INTERPOLATION_DEFAULT);
625};
626
627}}} // namespace plask::optical::effective
628
629#endif // PLASK__MODULE_OPTICAL_EFM_HPP