PLaSK library

Loading...
Searching...
No Matches
expansioncyl.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__SOLVER__SLAB_EXPANSIONCYL_H
15#define PLASK__SOLVER__SLAB_EXPANSIONCYL_H
16
17#include <plask/plask.hpp>
18
19#include "../expansion.hpp"
20#include "../meshadapter.hpp"
21#include "../patterson.hpp"
22
23 namespace plask { namespace optical { namespace modal {
24
25 struct BesselSolverCyl;
26
27 struct PLASK_SOLVER_API ExpansionBessel : public Expansion {
28 int m;
29
30 bool initialized;
31
32 bool m_changed;
33
36 OrderedAxis rbounds;
37
39 std::vector<double> kpts;
40
42 shared_ptr<RectangularMesh<2>> mesh;
43
48 ExpansionBessel(BesselSolverCyl* solver);
49
50 virtual ~ExpansionBessel() {}
51
53 void init1();
54
56 virtual void init2() = 0;
57
59 void init3();
60
62 virtual void reset();
63
64 bool diagonalQE(size_t l) const override { return diagonals[l]; }
65
66 size_t matrixSize() const override;
67
68 void prepareField() override;
69
70 void cleanupField() override;
71
72 LazyData<Vec<3, dcomplex>> getField(size_t layer,
73 const shared_ptr<const typename LevelsAdapter::Level>& level,
74 const cvector& E,
75 const cvector& H) override;
76
77 LazyData<Tensor3<dcomplex>> getMaterialEps(size_t layer,
78 const shared_ptr<const typename LevelsAdapter::Level>& level,
79 InterpolationMethod interp = INTERPOLATION_DEFAULT) override;
80
81 double integratePoyntingVert(const cvector& E, const cvector& H) override;
82
83 double integrateField(WhichField field, size_t layer, const cmatrix& TE, const cmatrix& TH,
84 const std::function<std::pair<dcomplex,dcomplex>(size_t, size_t)>& vertical) override;
85
86 private:
87 inline std::pair<double,double> getTC(size_t layer, size_t ri) {
88 double T = 0., W = 0., C = 0.;
89 for (size_t k = 0, v = ri * solver->verts->size(); k != mesh->vert()->size(); ++v, ++k) {
90 if (solver->stack[k] == layer) {
91 double w = (k == 0 || k == mesh->vert()->size() - 1) ? 1e-6 : solver->vbounds->at(k) - solver->vbounds->at(k - 1);
92 T += w * temperature[v];
93 C += w * carriers[v];
94 W += w;
95 }
96 }
97 T /= W;
98 C /= W;
99 return std::make_pair(T, C);
100 }
101
102 protected:
104 struct Segment {
105 double Z;
106 double D;
107 DataVector<double> weights;
108 };
109
111 std::vector<Segment> segments;
112
114 std::vector<bool> diagonals;
115
117 struct Integrals {
118 cmatrix V_k;
119 cmatrix TT;
120 cmatrix Tss;
121 cmatrix Tsp;
122 cmatrix Tps;
123 cmatrix Tpp;
124
125 void reset() {
126 V_k.reset();
127 TT.reset();
128 Tss.reset();
129 Tsp.reset();
130 Tps.reset();
131 Tpp.reset();
132 }
133
134 void reset(size_t N) {
135 V_k.reset(N, N);
136 TT.reset(2*N, 2*N);
137 size_t NN = N * N;
138 Tss.reset(N, N, TT.data());
139 Tsp.reset(N, N, TT.data() + NN);
140 Tps.reset(N, N, TT.data() + 2*NN);
141 Tpp.reset(N, N, TT.data() + 3*NN);
142 }
143 };
144
146 std::vector<Integrals> layers_integrals;
147
148 void beforeLayersIntegrals(dcomplex lam, dcomplex glam) override;
149
150 Tensor3<dcomplex> getEps(size_t layer, size_t ri, double r, double matz, double lam, double glam);
151
152 void layerIntegrals(size_t layer, double lam, double glam) override;
153
154 virtual void integrateParams(Integrals& integrals,
155 const dcomplex* datap, const dcomplex* datar, const dcomplex* dataz,
156 dcomplex datap0, dcomplex datar0, dcomplex dataz0) = 0;
157
158 virtual double fieldFactor(size_t i) = 0;
159
160 virtual cmatrix getHzMatrix(const cmatrix& Bz, cmatrix& Hz) = 0;
161
162 public:
163
164 void beforeGetEpsilon() override;
165
166 void afterGetEpsilon() override;
167
168 unsigned getM() const { return m; }
169 void setM(unsigned n) {
170 if (int(n) != m) {
171 write_debug("{0}: m changed from {1} to {2}", solver->getId(), m, n);
172 m = n;
173 solver->recompute_integrals = true;
174 solver->clearFields();
175 }
176 }
177
179 size_t idxs(size_t i) { return 2 * i; }
180
182 size_t idxp(size_t i) { return 2 * i + 1; }
183
184 #ifndef NDEBUG
185 cmatrix epsV_k(size_t layer);
186 cmatrix epsTss(size_t layer);
187 cmatrix epsTsp(size_t layer);
188 cmatrix epsTps(size_t layer);
189 cmatrix epsTpp(size_t layer);
190 #endif
191
193 virtual std::vector<double> getKpts() {
194 std::vector<double> res;
195 res.reserve(kpts.size());
196 double ib = 1. / rbounds[rbounds.size()-1];
197 for (double k: kpts) res.push_back(k * ib);
198 return res;
199 }
200};
201
202}}} // namespace plask::optical::modal
203
204#endif // PLASK__SOLVER__SLAB_EXPANSIONCYL_H