PLaSK library

Loading...
Searching...
No Matches
solver2d.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_OPTICAL_MODAL_FOURIER_REFLECTION_2D_H
15#define PLASK__SOLVER_OPTICAL_MODAL_FOURIER_REFLECTION_2D_H
16
17#include <plask/plask.hpp>
18
19#include "../solver.hpp"
20#include "expansion2d.hpp"
21
22 namespace plask { namespace optical { namespace modal {
23
27 struct PLASK_SOLVER_API FourierSolver2D: public ModalSolver<SolverWithMesh<Geometry2DCartesian,MeshAxis>> {
28
29 typedef ModalSolver<SolverWithMesh<Geometry2DCartesian,MeshAxis>> BaseType;
30
31 friend struct ExpansionPW2D;
32
33 std::string getClassName() const override { return "optical.Fourier2D"; }
34
36 enum What {
37 WHAT_WAVELENGTH,
38 WHAT_K0,
39 WHAT_NEFF,
40 WHAT_KTRAN,
41 WHAT_BETA
42 };
43
44
46 struct Mode {
47 Expansion::Component symmetry;
48 Expansion::Component polarization;
49 double lam0;
50 dcomplex k0;
51 dcomplex beta;
52 dcomplex ktran;
53 double power;
54 double tolx;
55
56 Mode(const ExpansionPW2D& expansion, double tolx):
57 symmetry(expansion.symmetry),
58 polarization(expansion.polarization),
59 lam0(expansion.lam0),
60 k0(expansion.k0),
61 beta(expansion.beta),
62 ktran(expansion.ktran),
63 power(1.),
64 tolx(tolx) {}
65
66 bool operator==(const Mode& other) const {
67 return is_equal(k0, other.k0) && is_equal(beta, other.beta) && is_equal(ktran, other.ktran)
68 && symmetry == other.symmetry && polarization == other.polarization &&
69 ((isnan(lam0) && isnan(other.lam0)) || lam0 == other.lam0)
70 ;
71 }
72
73 bool operator==(const ExpansionPW2D& other) const {
74 return is_equal(k0, other.k0) && is_equal(beta, other.beta) && is_equal(ktran, other.ktran)
75 && symmetry == other.symmetry && polarization == other.polarization &&
76 ((isnan(lam0) && isnan(other.lam0)) || lam0 == other.lam0)
77 ;
78 }
79
80 template <typename T>
81 bool operator!=(const T& other) const {
82 return !(*this == other);
83 }
84
85 private:
86
88 template <typename T>
89 bool is_equal(T a, T b) const {
90 return abs(a-b) <= tolx;
91 }
92 };
93
95 enum FourierType {
96 FOURIER_DISCRETE,
97 FOURIER_ANALYTIC
98 };
99
100 protected:
101
102 dcomplex beta,
103 ktran;
104
105 Expansion::Component symmetry;
106 Expansion::Component polarization;
107
109 size_t size;
110
111 void onInitialize() override;
112
113 void onInvalidate() override;
114
115 void computeIntegrals() override {
116 expansion.computeIntegrals();
117 }
118
120 int dct;
121
123 FourierType ftt;
124
125 public:
126
128 ExpansionPW2D expansion;
129
131 std::vector<Mode> modes;
132
133 void clearModes() override {
134 modes.clear();
135 }
136
137 bool setExpansionDefaults(bool with_k0=true) override {
138 bool changed = false;
139 if (expansion.getLam0() != getLam0()) { changed = true; expansion.setLam0(getLam0()); }
140 if (with_k0) {
141 if (expansion.getK0() != getK0()) { changed = true; expansion.setK0(getK0()); }
142 }
143 if (expansion.getBeta() != getBeta()) { changed = true; expansion.setBeta(getBeta()); }
144 if (expansion.getKtran() != getKtran()) { changed = true; expansion.setKtran(getKtran()); }
145 if (expansion.getSymmetry() != getSymmetry()) { changed = true; expansion.setSymmetry(getSymmetry()); }
146 if (expansion.getPolarization() != getPolarization()) { changed = true; expansion.setPolarization(getPolarization()); }
147 return changed;
148 }
149
151 size_t refine;
152
154 PML pml;
155
157 plask::optional<std::pair<double,double>> mirrors;
158
160 ProviderFor<ModeEffectiveIndex>::Delegate outNeff;
161
162 FourierSolver2D(const std::string& name="");
163
164 void loadConfiguration(XMLReader& reader, Manager& manager) override;
165
169 void setSimpleMesh() {
170 writelog(LOG_DETAIL, "Creating simple mesh");
171 setMesh(plask::make_shared<OrderedMesh1DSimpleGenerator>());
172 }
173
181 size_t findMode(FourierSolver2D::What what, dcomplex start);
182
184 size_t getSize() const { return size; }
186 void setSize(size_t n) {
187 size = n;
188 invalidate();
189 }
190
192 int getDCT() const { return dct; }
194 void setDCT(int n) {
195 if (n != 1 && n != 2)
196 throw BadInput(getId(), "bad DCT type (can be only 1 or 2)");
197 if (dct != n) {
198 dct = n;
199 if (symmetric() && ftt == FOURIER_DISCRETE) invalidate();
200 }
201 }
203 bool dct2() const { return dct == 2; }
204
205
207 FourierType getFourierType() const { return ftt; }
209 void setFourierType(FourierType ft) {
210 if (ft != ftt) {
211 ftt = ft;
212 invalidate();
213 }
214 }
215
217 dcomplex getKtran() const { return ktran; }
218
220 void setKtran(dcomplex k) {
221 if (k != 0. && (symmetric() || symmetry != Expansion::E_UNSPECIFIED)) {
222 Solver::writelog(LOG_WARNING, "Resetting mode symmetry");
223 symmetry = Expansion::E_UNSPECIFIED;
224 invalidate();
225 }
226 if (k != ktran && transfer) transfer->fields_determined = Transfer::DETERMINED_NOTHING;
227 ktran = k;
228 }
229
231 void setBeta(dcomplex k) {
232 if (k != 0. && (separated() || polarization != Expansion::E_UNSPECIFIED)) {
233 Solver::writelog(LOG_WARNING, "Resetting polarizations separation");
234 polarization = Expansion::E_UNSPECIFIED;
235 invalidate();
236 }
237 if (k != beta && transfer) transfer->fields_determined = Transfer::DETERMINED_NOTHING;
238 beta = k;
239 }
240
242 dcomplex getBeta() const { return beta; }
243
245 Expansion::Component getSymmetry() const { return symmetry; }
246
248 void setSymmetry(Expansion::Component sym) {
249 if (sym != Expansion::E_UNSPECIFIED && geometry && !geometry->isSymmetric(Geometry2DCartesian::DIRECTION_TRAN))
250 throw BadInput(getId(), "symmetry not allowed for asymmetric structure");
251 if ((symmetry == Expansion::E_UNSPECIFIED) != (sym == Expansion::E_UNSPECIFIED))
252 invalidate();
253 if (ktran != 0. && sym != Expansion::E_UNSPECIFIED) {
254 Solver::writelog(LOG_WARNING, "Resetting ktran to 0.");
255 ktran = 0.;
256 expansion.setKtran(0.);
257 }
258 symmetry = sym;
259 }
260
262 Expansion::Component getPolarization() const { return polarization; }
263
265 void setPolarization(Expansion::Component pol) {
266 if (polarization != pol)
267 invalidate();
268 if (beta != 0. && pol != Expansion::E_UNSPECIFIED) {
269 Solver::writelog(LOG_WARNING, "Resetting beta to 0.");
270 beta = 0.;
271 expansion.setBeta(0.);
272 }
273 polarization = pol;
274 }
275
277 bool symmetric() const { return symmetry != Expansion::E_UNSPECIFIED; }
278
280 bool separated() const { return polarization != Expansion::E_UNSPECIFIED; }
281
282 Expansion& getExpansion() override { return expansion; }
283
284 // RegularAxis getMesh() const { return *expansion.mesh->tran(); }
285
293 cvector incidentVector(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam=NAN);
294
295 using BaseType::incidentVector;
296
306 cvector incidentGaussian(Transfer::IncidentDirection side, Expansion::Component polarization,
307 double sigma, double center=0., dcomplex lam=NAN);
308
309 private:
310
311 size_t initIncidence(Transfer::IncidentDirection side, Expansion::Component polarization, dcomplex lam=NAN);
312
313 public:
314
323 LazyData<Vec<3,dcomplex>> getScatteredFieldE(const cvector& incident,
324 Transfer::IncidentDirection side,
325 shared_ptr<const MeshD<2>> dst_mesh,
326 InterpolationMethod method,
327 PropagationDirection part = PROPAGATION_TOTAL) {
328 if (!Solver::initCalculation()) setExpansionDefaults();
329 // expansion.setK0(2e3*PI / wavelength);
330 if (expansion.separated())
331 expansion.setPolarization(polarization);
332 if (!transfer) initTransfer(expansion, true);
333 return transfer->getScatteredFieldE(incident, side, dst_mesh, method, part);
334 }
335
344 LazyData<Vec<3,dcomplex>> getScatteredFieldH(const cvector& incident,
345 Transfer::IncidentDirection side,
346 shared_ptr<const MeshD<2>> dst_mesh,
347 InterpolationMethod method,
348 PropagationDirection part = PROPAGATION_TOTAL) {
349 if (!Solver::initCalculation()) setExpansionDefaults();
350 // expansion.setK0(2e3*PI / wavelength);
351 if (expansion.separated())
352 expansion.setPolarization(polarization);
353 if (!transfer) initTransfer(expansion, true);
354 return transfer->getScatteredFieldH(incident, side, dst_mesh, method, part);
355 }
356
364 LazyData<double> getScatteredFieldMagnitude(const cvector& incident,
365 Transfer::IncidentDirection side,
366 shared_ptr<const MeshD<2>> dst_mesh,
367 InterpolationMethod method) {
368 if (!Solver::initCalculation()) setExpansionDefaults();
369 // expansion.setK0(2e3*PI / wavelength);
370 if (expansion.separated())
371 expansion.setPolarization(polarization);
372 if (!transfer) initTransfer(expansion, true);
373 return transfer->getScatteredFieldMagnitude(incident, side, dst_mesh, method);
374 }
375
382 cvector getFieldVectorE(size_t num, double z) {
383 applyMode(modes[num]);
384 return transfer->getFieldVectorE(z);
385 }
386
393 cvector getFieldVectorH(size_t num, double z) {
394 applyMode(modes[num]);
395 return transfer->getFieldVectorH(z);
396 }
397
405 double getIntegralEE(size_t num, double z1, double z2) {
406 applyMode(modes[num]);
407 return transfer->getFieldIntegral(FIELD_E, z1, z2, modes[num].power);
408 }
409
417 double getIntegralHH(size_t num, double z1, double z2) {
418 applyMode(modes[num]);
419 return transfer->getFieldIntegral(FIELD_H, z1, z2, modes[num].power);
420 }
421
429 cvector getScatteredFieldVectorE(const cvector& incident, Transfer::IncidentDirection side, double z) {
430 if (!Solver::initCalculation()) setExpansionDefaults(false);
431 if (!transfer) initTransfer(expansion, true);
432 return transfer->getScatteredFieldVectorE(incident, side, z, PROPAGATION_TOTAL);
433 }
434
442 cvector getScatteredFieldVectorH(const cvector& incident, Transfer::IncidentDirection side, double z) {
443 if (!Solver::initCalculation()) setExpansionDefaults(false);
444 if (!transfer) initTransfer(expansion, true);
445 return transfer->getScatteredFieldVectorH(incident, side, z, PROPAGATION_TOTAL);
446 }
447
456 double getScatteredIntegralEE(const cvector& incident, Transfer::IncidentDirection side, double z1, double z2) {
457 if (!Solver::initCalculation()) setExpansionDefaults(false);
458 if (!transfer) initTransfer(expansion, true);
459 return transfer->getScatteredFieldIntegral(FIELD_E, incident, side, z1, z2);
460 }
461
470 double getScatteredIntegralHH(const cvector& incident, Transfer::IncidentDirection side, double z1, double z2) {
471 if (!Solver::initCalculation()) setExpansionDefaults(false);
472 if (!transfer) initTransfer(expansion, true);
473 return transfer->getScatteredFieldIntegral(FIELD_H, incident, side, z1, z2);
474 }
475
477 size_t setMode() {
478 if (abs2(this->getDeterminant()) > root.tolf_max*root.tolf_max)
479 throw BadInput(this->getId(), "cannot set the mode, determinant too large");
480 return insertMode();
481 }
482
483 protected:
484
486 size_t insertMode() {
487 static bool warn = true;
488 if (warn && emission != EMISSION_TOP && emission != EMISSION_BOTTOM) {
489 writelog(LOG_WARNING, "Mode fields are not normalized");
490 warn = false;
491 }
492 Mode mode(expansion, root.tolx);
493 for (size_t i = 0; i != modes.size(); ++i)
494 if (modes[i] == mode) return i;
495 modes.push_back(mode);
496 outNeff.fireChanged();
497 outLightMagnitude.fireChanged();
498 outLightE.fireChanged();
499 outLightH.fireChanged();
500 return modes.size()-1;
501 }
502
503 size_t nummodes() const override { return modes.size(); }
504
505 double applyMode(size_t n) override {
506 if (n >= modes.size()) throw BadInput(this->getId(), "mode {0} has not been computed", n);
507 applyMode(modes[n]);
508 return modes[n].power;
509 }
510
515 dcomplex getEffectiveIndex(size_t n) {
516 if (n >= modes.size()) throw NoValue(ModeEffectiveIndex::NAME);
517 return modes[n].beta / modes[n].k0;
518 }
519
521 double getMirrorLosses(double n) {
522 double L = geometry->getExtrusion()->getLength();
523 if (isinf(L)) return 0.;
524 const double lambda = real(2e3*PI / k0);
525 double R1, R2;
526 if (mirrors) {
527 std::tie(R1,R2) = *mirrors;
528 } else {
529 const double n1 = real(geometry->getFrontMaterial()->Nr(lambda, 300.)),
530 n2 = real(geometry->getBackMaterial()->Nr(lambda, 300.));
531 R1 = (n-n1) / (n+n1); R1 *= R1;
532 R2 = (n-n2) / (n+n2); R2 *= R2;
533 }
534 return 0.5 * std::log(R1*R2) / L;
535 }
536
537 void applyMode(const Mode& mode) {
538 writelog(LOG_DEBUG, "Current mode <lam: {:.2f}nm, neff: {}, ktran: {}/um, polarization: {}, symmetry: {}>",
539 real(2e3*PI/mode.k0),
540 str(mode.beta/mode.k0, "{:.3f}{:+.3g}j"),
541 str(mode.ktran, "({:.3g}{:+.3g}j)", "{:.3g}"),
542 (mode.polarization == Expansion::E_LONG)? "El" : (mode.polarization == Expansion::E_TRAN)? "Et" : "none",
543 (mode.symmetry == Expansion::E_LONG)? "El" : (mode.symmetry == Expansion::E_TRAN)? "Et" : "none"
544 );
545 if (mode != expansion) {
546 expansion.setLam0(mode.lam0);
547 expansion.setK0(mode.k0);
548 expansion.beta = mode.beta;
549 expansion.ktran = mode.ktran;
550 expansion.symmetry = mode.symmetry;
551 expansion.polarization = mode.polarization;
552 clearFields();
553 }
554 }
555
556 double getWavelength(size_t n) override;
557};
558
559
560}}} // namespace
561
562#endif