PLaSK library

Loading...
Searching...
No Matches
eim.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_EIM_HPP
15#define PLASK__MODULE_OPTICAL_EIM_HPP
16
17#include <limits>
18
19#include <plask/plask.hpp>
20
21#include "bisection.hpp"
22#include "rootdigger.hpp"
23
24 namespace plask { namespace optical { namespace effective {
25
29 struct PLASK_SOLVER_API EffectiveIndex2D : public SolverWithMesh<Geometry2DCartesian, RectangularMesh<2>> {
31 enum Symmetry { SYMMETRY_DEFAULT, SYMMETRY_POSITIVE, SYMMETRY_NEGATIVE, SYMMETRY_NONE };
32
34 enum Polarization {
35 TE,
36 TM,
37 };
38
40 enum Emission { FRONT, BACK };
41
42 struct Field {
43 dcomplex F, B;
44 Field() = default;
45 Field(dcomplex f, dcomplex b) : F(f), B(b) {}
46 Field operator*(dcomplex a) const { return Field(F * a, B * a); }
47 Field operator/(dcomplex a) const { return Field(F / a, B / a); }
48 Field operator*=(dcomplex a) {
49 F *= a;
50 B *= a;
51 return *this;
52 }
53 Field operator/=(dcomplex a) {
54 F /= a;
55 B /= a;
56 return *this;
57 }
58 };
59
60 struct Matrix {
61 dcomplex ff, fb, bf, bb;
62 Matrix() = default;
63 Matrix(dcomplex t1, dcomplex t2, dcomplex t3, dcomplex t4) : ff(t1), fb(t2), bf(t3), bb(t4) {}
64 static Matrix eye() { return Matrix(1., 0., 0., 1.); }
65 static Matrix diag(dcomplex f, dcomplex b) { return Matrix(f, 0., 0., b); }
66 Matrix operator*(const Matrix& T) {
67 return Matrix(ff * T.ff + fb * T.bf, ff * T.fb + fb * T.bb, bf * T.ff + bb * T.bf, bf * T.fb + bb * T.bb);
68 }
69 Field solve(const Field& v) { return Field(bb * v.F - fb * v.B, -bf * v.F + ff * v.B) / (ff * bb - fb * bf); }
70 };
71
73 struct Mode {
74 EffectiveIndex2D* solver;
75 Symmetry symmetry;
76 dcomplex neff;
77 bool have_fields;
78 std::vector<Field, aligned_allocator<Field>> xfields;
79 std::vector<double, aligned_allocator<double>> xweights;
80 double power;
81
82 Mode(EffectiveIndex2D* solver, Symmetry sym)
83 : solver(solver), have_fields(false), xfields(solver->xend), xweights(solver->xend), power(1.) {
84 setSymmetry(sym);
85 }
86
87 void setSymmetry(Symmetry sym) {
88 if (solver->geometry->isSymmetric(Geometry::DIRECTION_TRAN)) {
89 if (sym == SYMMETRY_DEFAULT)
90 sym = SYMMETRY_POSITIVE;
91 else if (sym == SYMMETRY_NONE)
92 throw BadInput(solver->getId(), "for symmetric geometry specify positive or negative symmetry");
93 } else {
94 if (sym == SYMMETRY_DEFAULT)
95 sym = SYMMETRY_NONE;
96 else if (sym != SYMMETRY_NONE)
97 throw BadInput(solver->getId(), "for non-symmetric geometry no symmetry may be specified");
98 }
99 symmetry = sym;
100 }
101
102 bool operator==(const Mode& other) const { return symmetry == other.symmetry && is_zero(neff - other.neff); }
103
105 double loss() const { return -2e7 * imag(neff * solver->k0); }
106 };
107
108 protected:
109 friend struct RootDigger;
110
111 size_t xbegin,
112 xend,
113 ybegin,
114 yend;
115
117 DataLog<dcomplex, dcomplex> log_value;
118
120 std::vector<std::vector<dcomplex, aligned_allocator<dcomplex>>> nrCache;
121
123 std::vector<Field, aligned_allocator<Field>> yfields;
124
126 std::vector<double, aligned_allocator<double>> yweights;
127
129 std::vector<dcomplex, aligned_allocator<dcomplex>> epsilons;
130
131 double stripex;
132
133 Polarization polarization;
134
135 bool recompute_neffs;
136
137 public:
138 Emission emission;
139
140 dcomplex vneff;
141
143 plask::optional<std::pair<double, double>> mirrors;
144
146 RootDigger::Params root;
147
149 RootDigger::Params stripe_root;
150
152 std::vector<Mode> modes;
153
155 ReceiverFor<Temperature, Geometry2DCartesian> inTemperature;
156
158 ReceiverFor<Gain, Geometry2DCartesian> inGain;
159
161 ReceiverFor<CarriersConcentration, Geometry2DCartesian> inCarriersConcentration;
162
164 typename ProviderFor<ModeEffectiveIndex>::Delegate outNeff;
165
167 typename ProviderFor<ModeLightMagnitude, Geometry2DCartesian>::Delegate outLightMagnitude;
168
170 typename ProviderFor<ModeLightE, Geometry2DCartesian>::Delegate outLightE;
171
173 typename ProviderFor<RefractiveIndex, Geometry2DCartesian>::Delegate outRefractiveIndex;
174
176 typename ProviderFor<Heat, Geometry2DCartesian>::Delegate outHeat;
177
178 EffectiveIndex2D(const std::string& name = "");
179
180 virtual ~EffectiveIndex2D() {
181 inTemperature.changedDisconnectMethod(this, &EffectiveIndex2D::onInputChange);
182 inGain.changedDisconnectMethod(this, &EffectiveIndex2D::onInputChange);
183 inCarriersConcentration.changedDisconnectMethod(this, &EffectiveIndex2D::onInputChange);
184 }
185
186 std::string getClassName() const override { return "optical.EffectiveIndex2D"; }
187
188 std::string getClassDescription() const override {
189 return "Calculate optical modes and optical field distribution using the effective index method "
190 "in Cartesian two-dimensional space.";
191 }
192
193 void loadConfiguration(plask::XMLReader& reader, plask::Manager& manager) override;
194
196 double getStripeX() const { return stripex; }
197
202 void setStripeX(double x) {
203 stripex = x;
204 invalidate();
205 }
206
208 Polarization getPolarization() const { return polarization; }
209
214 void setPolarization(Polarization polar) {
215 polarization = polar;
216 invalidate();
217 }
218
220 dcomplex getWavelength() const { return 2e3 * PI / k0; }
221
226 void setWavelength(dcomplex wavelength) {
227 k0 = 2e3 * PI / wavelength;
228 invalidate();
229 }
230
234 void setSimpleMesh() {
235 writelog(LOG_DETAIL, "Creating simple mesh");
236 setMesh(plask::make_shared<RectangularMesh2DSimpleGenerator>());
237 }
238
245 void setHorizontalMesh(shared_ptr<MeshAxis> meshx) { // TODO pointer to mesh is held now, is this fine?
246 writelog(LOG_DETAIL, "Setting horizontal mesh");
247 if (!geometry) throw NoChildException();
248 auto meshxy = RectangularMesh2DSimpleGenerator().generate_t<RectangularMesh<2>>(geometry->getChild());
249 meshxy->setAxis0(meshx);
250 setMesh(meshxy);
251 }
252
263 std::vector<dcomplex> searchVNeffs(plask::dcomplex neff1 = 0.,
264 plask::dcomplex neff2 = 0.,
265 size_t resteps = 256,
266 size_t imsteps = 64,
267 dcomplex eps = dcomplex(1e-6, 1e-9));
268
275 size_t findMode(dcomplex neff, Symmetry symmetry = SYMMETRY_DEFAULT);
276
287 std::vector<size_t> findModes(dcomplex neff1 = 0.,
288 dcomplex neff2 = 0.,
289 Symmetry symmetry = SYMMETRY_DEFAULT,
290 size_t resteps = 256,
291 size_t imsteps = 64,
292 dcomplex eps = dcomplex(1e-6, 1e-9));
293
298 dcomplex getVertDeterminant(dcomplex neff) {
299 updateCache();
300 size_t stripe = mesh->tran()->findIndex(stripex);
301 if (stripe < xbegin)
302 stripe = xbegin;
303 else if (stripe >= xend)
304 stripe = xend - 1;
305 return detS1(neff, nrCache[stripe]);
306 }
307
313 dcomplex getDeterminant(dcomplex neff, Symmetry sym = SYMMETRY_DEFAULT) {
314 stageOne();
315 Mode mode(this, sym);
316 dcomplex det = detS(neff, mode);
317 return det;
318 }
319
327 size_t setMode(dcomplex neff, Symmetry sym = SYMMETRY_DEFAULT);
328
330 void clearModes() { modes.clear(); }
331
336 double getTotalAbsorption(Mode& mode);
337
342 double getTotalAbsorption(size_t num);
343
348 dcomplex getDeltaNeff(double x) {
349 stageOne();
350 size_t i(clamp(mesh->tran()->findIndex(geometry->isSymmetric(Geometry::DIRECTION_TRAN) ? abs(x) : x), xbegin, xend - 1));
351 return sqrt(epsilons[i]);
352 }
353
354 protected:
356 void onInputChange(ReceiverBase&, ReceiverBase::ChangeReason) { invalidate(); }
357
359 void onInitialize() override;
360
362 void onInvalidate() override;
363
365 dcomplex k0;
366
368 bool need_gain;
369
371 double getMirrorLosses(dcomplex n) {
372 double L = geometry->getExtrusion()->getLength();
373 if (isinf(L)) return 0.;
374 const double lambda = real(2e3 * PI / k0);
375 double R1, R2;
376 if (mirrors) {
377 std::tie(R1, R2) = *mirrors;
378 } else {
379 const double n1 = real(geometry->getFrontMaterial()->Nr(lambda, 300.)),
380 n2 = real(geometry->getBackMaterial()->Nr(lambda, 300.));
381 R1 = abs((n - n1) / (n + n1));
382 R2 = abs((n - n2) / (n + n2));
383 }
384 return lambda * std::log(R1 * R2) / (4e3 * PI * L);
385 }
386
391 void updateCache();
392
397 void stageOne();
398
403 void computeWeights(size_t stripe);
404
409 void normalizeFields(Mode& mode, const std::vector<dcomplex, aligned_allocator<dcomplex>>& kx);
410
417 dcomplex detS1(const dcomplex& x, const std::vector<dcomplex, aligned_allocator<dcomplex>>& NR, bool save = false);
418
424 dcomplex detS(const dcomplex& x, Mode& mode, bool save = false);
425
427 size_t insertMode(const Mode& mode) {
428 for (size_t i = 0; i != modes.size(); ++i)
429 if (modes[i] == mode) return i;
430 modes.push_back(mode);
431 outNeff.fireChanged();
432 outLightMagnitude.fireChanged();
433 outLightE.fireChanged();
434 return modes.size() - 1;
435 }
436
438 size_t nmodes() const { return modes.size(); }
439
444 dcomplex getEffectiveIndex(size_t n) {
445 if (n >= modes.size()) throw NoValue(ModeEffectiveIndex::NAME);
446 return modes[n].neff;
447 }
448
449 template <typename T> struct FieldDataBase;
450 template <typename T> struct FieldDataInefficient;
451 template <typename T> struct FieldDataEfficient;
452 struct HeatDataImpl;
453
455 const LazyData<double> getLightMagnitude(std::size_t num,
456 shared_ptr<const plask::MeshD<2>> dst_mesh,
457 plask::InterpolationMethod = INTERPOLATION_DEFAULT);
458
460 const LazyData<Vec<3, dcomplex>> getElectricField(std::size_t num,
461 shared_ptr<const plask::MeshD<2>> dst_mesh,
462 plask::InterpolationMethod = INTERPOLATION_DEFAULT);
463
465 const LazyData<dcomplex> getRefractiveIndex(RefractiveIndex::EnumType component,
466 shared_ptr<const MeshD<2>> dst_mesh,
467 dcomplex lam,
468 InterpolationMethod = INTERPOLATION_DEFAULT);
469
471 const LazyData<double> getHeat(shared_ptr<const MeshD<2>> dst_mesh, InterpolationMethod method = INTERPOLATION_DEFAULT);
472};
473
474}}} // namespace plask::optical::effective
475
476#endif // PLASK__MODULE_OPTICAL_EIM_HPP