PLaSK library

Loading...
Searching...
No Matches
solvercyl.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_SOLVERCYL_H
15#define PLASK__SOLVER__SLAB_SOLVERCYL_H
16
17#include <plask/plask.hpp>
18
19#include "../solver.hpp"
20#include "../reflection.hpp"
21#include "expansioncyl-fini.hpp"
22#include "expansioncyl-infini.hpp"
23
24
25 namespace plask { namespace optical { namespace modal {
26
30 struct PLASK_SOLVER_API BesselSolverCyl: public ModalSolver<SolverWithMesh<Geometry2DCylindrical,MeshAxis>> {
31
32 typedef ModalSolver<SolverWithMesh<Geometry2DCylindrical,MeshAxis>> BaseType;
33
34 friend struct ExpansionBessel;
35 friend struct ExpansionBesselFini;
36 friend struct ExpansionBesselInfini;
37
38 enum BesselDomain {
39 DOMAIN_FINITE,
40 DOMAIN_INFINITE
41 };
42
43 enum InfiniteWavevectors {
44 WAVEVECTORS_UNIFORM,
45 WAVEVECTORS_NONUNIFORM,
46 // WAVEVECTORS_LEGENDRE,
47 WAVEVECTORS_LAGUERRE,
48 WAVEVECTORS_MANUAL
49 };
50
51 enum Rule {
52 RULE_DIRECT,
53 RULE_COMBINED_1,
54 RULE_COMBINED_2,
55 RULE_OLD
56 };
57
58 const char* ruleName() {
59 switch (rule) {
60 case RULE_DIRECT: return "inverse";
61 case RULE_COMBINED_1: return "inverse (mod 1)";
62 case RULE_COMBINED_2: return "inverse (mod 2)";
63 case RULE_OLD: return "direct";
64 }
65 return "unknown";
66 }
67
68 std::string getClassName() const override { return "optical.BesselCyl"; }
69
70 struct Mode {
71 double lam0;
72 dcomplex k0;
73 int m;
74 double power;
75 double tolx;
76
77 Mode(const std::unique_ptr<ExpansionBessel>& expansion, double tolx):
78 lam0(expansion->lam0), k0(expansion->k0), m(expansion->m), power(1.), tolx(tolx) {}
79
80 bool operator==(const Mode& other) const {
81 return m == other.m && is_equal(k0, other.k0) && is_equal(lam0, other.lam0) &&
82 ((isnan(lam0) && isnan(other.lam0)) || lam0 == other.lam0);
83 }
84
85 bool operator==(const std::unique_ptr<ExpansionBessel>& other) const {
86 return m == other->m && is_equal(k0, other->k0) && is_equal(lam0, other->lam0) &&
87 ((isnan(lam0) && isnan(other->lam0)) || lam0 == other->lam0);
88 }
89
90 template <typename T>
91 bool operator!=(const T& other) const {
92 return !(*this == other);
93 }
94
95 private:
96
98 template <typename T>
99 bool is_equal(T a, T b) const {
100 return abs(a-b) <= tolx;
101 }
102 };
103
104 protected:
105
107 BesselDomain domain;
108
110 int m;
111
113 size_t size;
114
115 void onInitialize() override;
116
117 void onInvalidate() override;
118
119 void computeIntegrals() override {
120 expansion->computeIntegrals();
121 }
122
124 Rule rule;
125
127 double kscale;
128
130 double kmax;
131
133 InfiniteWavevectors kmethod;
134
135 public:
136
138 std::vector<double> klist;
139
141 boost::optional<std::vector<double>> kweights;
142
144 std::unique_ptr<ExpansionBessel> expansion;
145
147 std::vector<Mode> modes;
148
150 std::vector<double> getKlist() {
151 initCalculation();
152 computeIntegrals();
153 return expansion->getKpts();
154 }
155
157 void setKlist(const std::vector<double>& values) {
158 if (kmethod != WAVEVECTORS_MANUAL) {
159 invalidate();
160 this->writelog(LOG_WARNING, "Setting Hankel transform method to Manual");
161 kmethod = WAVEVECTORS_MANUAL;
162 }
163 klist = values;
164 }
165
167 std::vector<double> getKweights() {
168 if (domain == DOMAIN_INFINITE) {
169 initCalculation();
170 computeIntegrals();
171 ExpansionBesselInfini* ex = dynamic_cast<ExpansionBesselInfini*>(expansion.get());
172 if (ex) return std::vector<double>(ex->kdelts.begin(), ex->kdelts.end());
173 }
174 return std::vector<double>();
175 }
176
178 void setKweights(const std::vector<double>& values) {
179 if (kmethod != WAVEVECTORS_MANUAL) {
180 invalidate();
181 this->writelog(LOG_WARNING, "Setting Hankel transform method to Manual");
182 kmethod = WAVEVECTORS_MANUAL;
183 }
184 kweights.reset(values);
185 }
186
187 void clearKweights() {
188 kweights.reset();
189 }
190
191 void clearModes() override {
192 modes.clear();
193 }
194
195 bool setExpansionDefaults(bool with_k0=true) override {
196 bool changed = false;
197 if (expansion->getLam0() != getLam0()) { changed = true; expansion->setLam0(getLam0()); }
198 if (with_k0) {
199 if (expansion->getK0() != getK0()) { changed = true; expansion->setK0(getK0()); }
200 }
201 if (expansion->getM() != getM()) { changed = true; expansion->setM(getM()); }
202 return changed;
203 }
204
206 double integral_error;
207
209 size_t max_integration_points;
210
212 PML pml;
213
215 typename ProviderFor<ModeLoss>::Delegate outLoss;
216
217 BesselSolverCyl(const std::string& name="");
218
219 void loadConfiguration(XMLReader& reader, Manager& manager) override;
220
227 size_t findMode(dcomplex start, int m=1);
228
230 size_t getSize() const { return size; }
232 void setSize(size_t n) {
233 size = n;
234 invalidate();
235 }
236
238 double getKscale() const { return kscale; }
240 void setKscale(double s) { kscale = s; }
241
243 double getKmax() const { return kmax; }
245 void setKmax(double s) { kmax = s; }
246
248 InfiniteWavevectors getKmethod() const { return kmethod; }
250 void setKmethod(InfiniteWavevectors k) { kmethod = k; }
251
253 BesselDomain getDomain() const { return domain; }
255 void setDomain(BesselDomain dom) {
256 domain = dom;
257 invalidate();
258 }
259
261 Rule getRule() const { return rule; }
263 void setRule(Rule r) {
264 rule = r;
265 invalidate();
266 }
267
269 unsigned getM() const { return m; }
271 void setM(unsigned n) { m = n; }
272
273 Expansion& getExpansion() override { return *expansion; }
274
283 LazyData<Vec<3,dcomplex>> getScatteredFieldE(const cvector& incident,
284 Transfer::IncidentDirection side,
285 const shared_ptr<const MeshD<2>>& dst_mesh,
286 InterpolationMethod method,
287 PropagationDirection part = PROPAGATION_TOTAL) {
288 if (!Solver::initCalculation()) setExpansionDefaults(false);
289 if (!transfer) initTransfer(*expansion, true);
290 return transfer->getScatteredFieldE(incident, side, dst_mesh, method, part);
291 }
292
301 LazyData<Vec<3,dcomplex>> getScatteredFieldH(const cvector& incident,
302 Transfer::IncidentDirection side,
303 const shared_ptr<const MeshD<2>>& dst_mesh,
304 InterpolationMethod method,
305 PropagationDirection part = PROPAGATION_TOTAL) {
306 if (!Solver::initCalculation()) setExpansionDefaults(false);
307 if (!transfer) initTransfer(*expansion, true);
308 return transfer->getScatteredFieldH(incident, side, dst_mesh, method, part);
309 }
310
318 LazyData<double> getScatteredFieldMagnitude(const cvector& incident,
319 Transfer::IncidentDirection side,
320 const shared_ptr<const MeshD<2>>& dst_mesh,
321 InterpolationMethod method) {
322 if (!Solver::initCalculation()) setExpansionDefaults(false);
323 if (!transfer) initTransfer(*expansion, true);
324 return transfer->getScatteredFieldMagnitude(incident, side, dst_mesh, method);
325 }
326
333 cvector getFieldVectorE(size_t num, double z) {
334 applyMode(modes[num]);
335 return transfer->getFieldVectorE(z);
336 }
337
344 cvector getFieldVectorH(size_t num, double z) {
345 applyMode(modes[num]);
346 return transfer->getFieldVectorH(z);
347 }
348
356 double getIntegralEE(size_t num, double z1, double z2) {
357 applyMode(modes[num]);
358 return transfer->getFieldIntegral(FIELD_E, z1, z2, modes[num].power);
359 }
360
368 double getIntegralHH(size_t num, double z1, double z2) {
369 applyMode(modes[num]);
370 return transfer->getFieldIntegral(FIELD_H, z1, z2, modes[num].power);
371 }
372
380 cvector getScatteredFieldVectorE(const cvector& incident, Transfer::IncidentDirection side, double z) {
381 if (!Solver::initCalculation()) setExpansionDefaults(false);
382 if (!transfer) initTransfer(*expansion, true);
383 return transfer->getScatteredFieldVectorE(incident, side, z, PROPAGATION_TOTAL);
384 }
385
393 cvector getScatteredFieldVectorH(const cvector& incident, Transfer::IncidentDirection side, double z) {
394 if (!Solver::initCalculation()) setExpansionDefaults(false);
395 if (!transfer) initTransfer(*expansion, true);
396 return transfer->getScatteredFieldVectorH(incident, side, z, PROPAGATION_TOTAL);
397 }
398
407 double getScatteredIntegralEE(const cvector& incident, Transfer::IncidentDirection side, double z1, double z2) {
408 if (!Solver::initCalculation()) setExpansionDefaults(false);
409 if (!transfer) initTransfer(*expansion, true);
410 return transfer->getScatteredFieldIntegral(FIELD_E, incident, side, z1, z2);
411 }
412
421 double getScatteredIntegralHH(const cvector& incident, Transfer::IncidentDirection side, double z1, double z2) {
422 if (!Solver::initCalculation()) setExpansionDefaults(false);
423 if (!transfer) initTransfer(*expansion, true);
424 return transfer->getScatteredFieldIntegral(FIELD_H, incident, side, z1, z2);
425 }
426
428 size_t setMode() {
429 if (abs2(this->getDeterminant()) > root.tolf_max*root.tolf_max)
430 throw BadInput(this->getId(), "cannot set the mode, determinant too large");
431 return insertMode();
432 }
433
434 protected:
436 size_t insertMode() {
437 static bool warn = true;
438 if (warn && ((emission != EMISSION_TOP && emission != EMISSION_BOTTOM) || domain == DOMAIN_INFINITE)) {
439 if (domain == DOMAIN_INFINITE)
440 writelog(LOG_WARNING, "Mode fields are not normalized (infinite domain)");
441 else
442 writelog(LOG_WARNING, "Mode fields are not normalized (emission direction not specified)");
443 warn = false;
444 }
445 Mode mode(expansion, root.tolx);
446 for (size_t i = 0; i != modes.size(); ++i)
447 if (modes[i] == mode) return i;
448 modes.push_back(mode);
449 outWavelength.fireChanged();
450 outLoss.fireChanged();
451 outLightMagnitude.fireChanged();
452 outLightE.fireChanged();
453 outLightH.fireChanged();
454 return modes.size()-1;
455 }
456
457 size_t nummodes() const override { return modes.size(); }
458
459 double applyMode(size_t n) override {
460 if (n >= modes.size()) throw BadInput(this->getId(), "mode {0} has not been computed", n);
461 applyMode(modes[n]);
462 return modes[n].power;
463 }
464
465 void applyMode(const Mode& mode) {
466 writelog(LOG_DEBUG, "Current mode <m: {:d}, lam: {}nm>", mode.m, str(2e3*PI/mode.k0, "({:.3f}{:+.3g}j)"));
467 expansion->setLam0(mode.lam0);
468 expansion->setK0(mode.k0);
469 expansion->setM(mode.m);
470 }
471
476 double getModalLoss(size_t n) {
477 if (n >= modes.size()) throw NoValue(ModeLoss::NAME);
478 return 2e4 * modes[n].k0.imag(); // 2e4 2/µm -> 2/cm
479 }
480
481 double getWavelength(size_t n) override;
482
483#ifndef NDEBUG
484 public:
485 cmatrix epsV_k(size_t layer);
486 cmatrix epsTss(size_t layer);
487 cmatrix epsTsp(size_t layer);
488 cmatrix epsTps(size_t layer);
489 cmatrix epsTpp(size_t layer);
490 cmatrix muV_k();
491 cmatrix muTss();
492 cmatrix muTsp();
493 cmatrix muTps();
494 cmatrix muTpp();
495#endif
496
497};
498
499
500
501}}} // # namespace plask::optical::modal
502
503#endif // PLASK__SOLVER__SLAB_SOLVERCYL_H