PLaSK library
Loading...
Searching...
No Matches
xance.cpp
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/*#if defined(_WIN32) || defined(__WIN32__) || defined(WIN32)
15# define BOOST_USE_WINDOWS_H
16#endif*/
17
18#include "xance.hpp"
19#include "solver.hpp"
20#include "diagonalizer.hpp"
21#include "expansion.hpp"
22#include "fortran.hpp"
23#include "meshadapter.hpp"
24
25namespace plask { namespace optical { namespace modal {
26
27XanceTransfer::XanceTransfer(ModalBase* solver, Expansion& expansion): Transfer(solver, expansion)
28{
29 // Reserve space for matrix multiplications...
30 const std::size_t N = diagonalizer->matrixSize();
31 Y = cmatrix(N,N);
32 needAllY = false;
33}
34
35
36void XanceTransfer::storeY(size_t n)
37{
38 if (needAllY) {
39 const std::size_t N = diagonalizer->matrixSize();
40 if (memY.size() != solver->stack.size()) {
41 // Allocate the storage for admittance matrices
42 memY.resize(solver->stack.size());
43 for (std::size_t i = 0; i < solver->stack.size(); i++) memY[i] = cmatrix(N,N);
44 }
45 memcpy(memY[n].data(), Y.data(), N*N*sizeof(dcomplex));
46 }
47}
48
49
50cvector XanceTransfer::getFieldVectorE(double z, std::size_t n, PropagationDirection part)
51{
52 cvector E0 = fields[n].E0;
53 cvector Ed = fields[n].Ed;
54
55 cdiagonal gamma = diagonalizer->Gamma(solver->stack[n]);
56 double d = get_d(n, z, part);
57
58 if ((n == 0 || std::size_t(n) == solver->vbounds->size()) && z < 0.)
59 return cvector(diagonalizer->source()->matrixSize(), NAN);
60
61 const std::size_t N = gamma.size();
62 cvector E(N);
63
64 switch (part) {
65 case PROPAGATION_TOTAL:
66 for (std::size_t i = 0; i < N; i++) {
67 dcomplex g = gamma[i];
68 //E[i] = (sin(g*(d-z)) * E0[i] + sin(g*z) * Ed[i]) / sin(g*d);
69 double a = abs(exp(2.*I*g*d));
70 if (isinf(a) || a < SMALL) {
71 dcomplex d0p = exp(I*g*z) - exp(I*g*(z-2*d));
72 dcomplex d0n = exp(I*g*(2*d-z)) - exp(-I*g*z);
73 if (isinf(real(d0p)) || isinf(imag(d0p))) d0p = 0.; else d0p = 1./ d0p;
74 if (isinf(real(d0n)) || isinf(imag(d0n))) d0n = 0.; else d0n = 1./ d0n;
75 dcomplex ddp = exp(I*g*(d-z)) - exp(-I*g*(d+z));
76 dcomplex ddn = exp(I*g*(d+z)) - exp(I*g*(z-d));
77 if (isinf(real(ddp)) || isinf(imag(ddp))) ddp = 0.; else ddp = 1./ ddp;
78 if (isinf(real(ddn)) || isinf(imag(ddn))) ddn = 0.; else ddn = 1./ ddn;
79 E[i] = (d0p-d0n) * E0[i] + (ddp-ddn) * Ed[i];
80 } else {
81 E[i] = (sinh(I*g*(d-z)) * E0[i] + sinh(I*g*z) * Ed[i]) / sinh(I*g*d);
82 }
83 }
84 break;
85 case PROPAGATION_UPWARDS:
86 for (std::size_t i = 0; i < N; i++) {
87 dcomplex g = gamma[i]; if (g.real() < 0) g = -g;
88 E[i] = 0.5 * (E0[i] * exp(I*g*d) - Ed[i]) * exp(-I*g*z) / sinh(I*g*d);
89 }
90 break;
91 case PROPAGATION_DOWNWARDS:
92 for (std::size_t i = 0; i < N; i++) {
93 dcomplex g = gamma[i]; if (g.real() < 0) g = -g;
94 E[i] = 0.5 * (Ed[i] - E0[i] * exp(-I*g*d)) * exp(I*g*z) / sinh(I*g*d);
95 }
96 break;
97 }
98
99 cvector result(diagonalizer->source()->matrixSize());
100 // result = diagonalizer->TE(n) * E;
101 mult_matrix_by_vector(diagonalizer->TE(solver->stack[n]), E, result);
102 return result;
103}
104
105
106cvector XanceTransfer::getFieldVectorH(double z, std::size_t n, PropagationDirection part)
107{
108 cvector H0 = fields[n].H0;
109 cvector Hd = fields[n].Hd;
110
111 cdiagonal gamma = diagonalizer->Gamma(solver->stack[n]);
112 double d = get_d(n, z, part);
113
114 if ((n == 0 || std::size_t(n) == solver->vbounds->size()) && z < 0.)
115 return cvector(diagonalizer->source()->matrixSize(), NAN);
116
117 const std::size_t N = gamma.size();
118 cvector H(N);
119
120 switch (part) {
121 case PROPAGATION_TOTAL:
122 for (std::size_t i = 0; i < N; i++) {
123 dcomplex g = gamma[i];
124 //H[i] = (sin(g*(d-z)) * H0[i] + sin(g*z) * Hd[i]) / sin(g*d);
125
126 double a = abs(exp(2.*I*g*d));
127 if (isinf(a) || a < SMALL) {
128 dcomplex d0p = exp(I*g*z) - exp(I*g*(z-2*d));
129 dcomplex d0n = exp(I*g*(2*d-z)) - exp(-I*g*z);
130 if (isinf(real(d0p)) || isinf(imag(d0p))) d0p = 0.; else d0p = 1./ d0p;
131 if (isinf(real(d0n)) || isinf(imag(d0n))) d0n = 0.; else d0n = 1./ d0n;
132 dcomplex ddp = exp(I*g*(d-z)) - exp(-I*g*(d+z));
133 dcomplex ddn = exp(I*g*(d+z)) - exp(I*g*(z-d));
134 if (isinf(real(ddp)) || isinf(imag(ddp))) ddp = 0.; else ddp = 1./ ddp;
135 if (isinf(real(ddn)) || isinf(imag(ddn))) ddn = 0.; else ddn = 1./ ddn;
136 H[i] = (d0p-d0n) * H0[i] + (ddp-ddn) * Hd[i];
137 } else {
138 H[i] = (sinh(I*g*(d-z)) * H0[i] + sinh(I*g*z) * Hd[i]) / sinh(I*g*d);
139 }
140 }
141 break;
142 case PROPAGATION_UPWARDS:
143 for (std::size_t i = 0; i < N; i++) {
144 dcomplex g = gamma[i]; if (g.real() < 0) g = -g;
145 H[i] = 0.5 * (H0[i] * exp(I*g*d) - Hd[i]) * exp(-I*g*z) / sinh(I*g*d);
146 }
147 break;
148 case PROPAGATION_DOWNWARDS:
149 for (std::size_t i = 0; i < N; i++) {
150 dcomplex g = gamma[i]; if (g.real() < 0) g = -g;
151 H[i] = 0.5 * (Hd[i] - H0[i] * exp(-I*g*d)) * exp(I*g*z) / sinh(I*g*d);
152 }
153 break;
154 }
155
156 cvector result(diagonalizer->source()->matrixSize());
157 // result = diagonalizer->TH(n) * H;
158 mult_matrix_by_vector(diagonalizer->TH(solver->stack[n]), H, result);
159 return result;
160}
161
162
163cvector XanceTransfer::getTransmissionVector(const cvector& incident, IncidentDirection side)
164{
165 determineReflectedFields(incident, side);
166 size_t n = (side == INCIDENCE_BOTTOM)? solver->stack.size()-1 : 0;
167 return fields[n].E0;
168}
169
170
171inline static dcomplex _int_xance(double z1, double z2, dcomplex a, dcomplex b = 0.) {
172 if (is_zero(a)) return 0.5 * (z2 - z1) * cosh(b);
173 const dcomplex a2 = 0.5 * a;
174 return cosh(a2 * (z1 + z2) + b) * sinh(a2 * (z2 - z1)) / a;
175}
176
177double XanceTransfer::integrateField(WhichField field, size_t n, double z1, double z2) {
178 size_t layer = solver->stack[n];
179 size_t N = diagonalizer->matrixSize();
180
181 const cvector& E0 = fields[n].E0;
182 const cvector& Ed = fields[n].Ed;
183 const cvector& H0 = fields[n].H0;
184 const cvector& Hd = fields[n].Hd;
185
186 cmatrix TE = diagonalizer->TE(layer),
187 TH = diagonalizer->TH(layer);
188 cdiagonal gamma = diagonalizer->Gamma(layer);
189
190 double d = get_d(n, z1, z2);
191
192 return diagonalizer->source()->integrateField(field, layer, TE, TH,
193 [z1, z2, d, gamma, E0, Ed, H0, Hd](size_t i, size_t j) -> std::pair<dcomplex,dcomplex> {
194 const dcomplex igm = I * (gamma[i] - conj(gamma[j])), igp = I * (gamma[i] + conj(gamma[j]));
195 const dcomplex igid = I * gamma[i] * d, igjd = I * conj(gamma[j]) * d;
196 dcomplex E = 0.;
197 dcomplex H = 0.;
198 if (!((is_zero(E0[i]) || is_zero(E0[j])) && (is_zero(H0[i]) || is_zero(H0[j])))) {
199 dcomplex val = _int_xance(z1, z2, -igm, igm*d) - _int_xance(z1, z2, -igp, igp*d);
200 E += E0[i] * conj(E0[j]) * val; H += H0[i] * conj(H0[j]) * val;
201 }
202 if (!((is_zero(E0[i]) || is_zero(Ed[j])) && (is_zero(H0[i]) || is_zero(Hd[j])))) {
203 dcomplex val = _int_xance(z1, z2, -igp, igid) - _int_xance(z1, z2, -igm, igid);
204 E += E0[i] * conj(Ed[j]) * val; H += H0[i] * conj(Hd[j]) * val;
205 }
206 if (!((is_zero(Ed[i]) || is_zero(E0[j])) && (is_zero(Hd[i]) || is_zero(H0[j])))) {
207 dcomplex val = _int_xance(z1, z2, igp, -igjd) - _int_xance(z1, z2, igm, igjd);
208 E += Ed[i] * conj(E0[j]) * val; H += Hd[i] * conj(H0[j]) * val;
209 }
210 if (!((is_zero(Ed[i]) || is_zero(Ed[j])) && (is_zero(Hd[i]) || is_zero(Hd[j])))) {
211 dcomplex val = _int_xance(z1, z2, igm) - _int_xance(z1, z2, igp);
212 E += Ed[i] * conj(Ed[j]) * val; H += Hd[i] * conj(Hd[j]) * val;
213 }
214 dcomplex f = 1. / (sinh(igid) * sinh(-igjd));
215 return std::make_pair(f * E, f * H);
216 });
217}
218
219}}} // namespace plask::optical::modal