PLaSK library

Loading...
Searching...
No Matches
reflection.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 "reflection.hpp"
19#include "solver.hpp"
20#include "diagonalizer.hpp"
21#include "expansion.hpp"
22#include "fortran.hpp"
23#include "meshadapter.hpp"
24
25#include <plask/utils/openmp.hpp>
26
27 namespace plask { namespace optical { namespace modal {
28
29 ReflectionTransfer::ReflectionTransfer(ModalBase* solver, Expansion& expansion, Matching matching): Transfer(solver, expansion),
30 matching(matching),
31 storeP(STORE_NONE) {
32 writelog(LOG_DETAIL, "{}: Initializing Reflection Transfer (with {} matching)", solver->getId(),
33 (matching == MATCH_ADMITTANCE)? "admittance" : "impedance");
34 size_t N = diagonalizer->matrixSize();
35 P = cmatrix(N,N);
36 phas = cdiagonal(N);
37 ipiv = aligned_new_array<int>(N);
38}
39
40
41 ReflectionTransfer::~ReflectionTransfer() {
42 // here we just use `aligned_free` as int is a PODT and diagonalizer->matrixSize() will call pure virtual method
43 // size_t N = diagonalizer->matrixSize();
44 // aligned_delete_array<int>(N, ipiv); ipiv = nullptr;
45 aligned_free(ipiv);
46}
47
48
49 void ReflectionTransfer::getFinalMatrix() {
50 getAM(0, solver->interface-1, false);
51 getAM(solver->stack.size()-1, solver->interface, true);
52}
53
54
55 void ReflectionTransfer::getAM(size_t start, size_t end, bool add, double mfac)
56{
57 // Get matrices sizes
58 const std::size_t N0 = diagonalizer->source()->matrixSize();
59 const std::size_t N = diagonalizer->matrixSize(); // <= N0
60 const std::size_t NN = N*N;
61 cmatrix work(N, N0, wrk); // matrix object for the workspace
62
63 findReflection(start, end, false, int(add)&1);
64
65 cdiagonal gamma = diagonalizer->Gamma(solver->stack[end]);
66
67 double H = (end == 0 || end == solver->vbounds->size())?
68 0 : abs(solver->vbounds->at(end) - solver->vbounds->at(end-1));
69 for (std::size_t i = 0; i < N; i++) phas[i] = exp(-I*gamma[i]*H);
70
71 mult_diagonal_by_matrix(phas, P); mult_matrix_by_diagonal(P, phas); // P = phas * P * phas
72 memcpy(temp.data(), P.data(), NN*sizeof(dcomplex)); // temp = P
73
74 double II = (matching == MATCH_ADMITTANCE)? 1. : -1.;
75
76 // temp = [ phas*P*phas ∓ I ] [ phas*P*phas ± I ]^{-1}
77 for (std::size_t i = 0, ii = 0; i < N; i++, ii += (N+1)) P[ii] += II; // P = P ± I
78 for (std::size_t i = 0, ii = 0; i < N; i++, ii += (N+1)) temp[ii] -= II; // temp = temp ∓ I
79 int info;
80 zgetrf(int(N), int(N), P.data(), int(N), ipiv, info); // P = LU(P)
81 ztrsm('R', 'U', 'N', 'N', int(N), int(N), 1., P.data(), int(N), temp.data(), int(N)); // temp = temp * U^{-1}
82 ztrsm('R', 'L', 'N', 'U', int(N), int(N), 1., P.data(), int(N), temp.data(), int(N)); // temp = temp * L^{-1}
83 // reorder columns (there is no such function in LAPACK)
84 for (std::ptrdiff_t j = N-1; j >=0; j--) {
85 int jp = ipiv[j]-1;
86 for (std::size_t i = 0; i < N; i++) std::swap(temp(i,j), temp(i,jp));
87 }
88
89 // M for the half of the structure
90 if (matching == MATCH_ADMITTANCE) {
91 mult_matrix_by_matrix(temp, diagonalizer->invTE(solver->stack[end]), work); // work = temp * invTE[end]
92 zgemm('N','N', int(N0), int(N0), int(N), mfac, diagonalizer->TH(solver->stack[end]).data(), int(N0),
93 wrk, int(N), add?1.:0., M.data(), int(N0)); // M = mfac * TH[end] * work
94 } else {
95 mult_matrix_by_matrix(temp, diagonalizer->invTH(solver->stack[end]), work); // work = temp * invTH[end]
96 zgemm('N','N', int(N0), int(N0), int(N), mfac, diagonalizer->TE(solver->stack[end]).data(), int(N0),
97 wrk, int(N), add?1.:0., M.data(), int(N0)); // M = mfac * TE[end] * work
98 }
99}
100
101
102 void ReflectionTransfer::findReflection(std::size_t start, std::size_t end, bool emitting, int store)
103{
104 // Should be called from 0 to interface-1
105 // and from count-1 to interface
106
107 write_debug("{}: searching for reflection for layers {:d} to {:d}", solver->getId(), start, end);
108
109 const std::ptrdiff_t inc = (start < end) ? 1 : -1;
110
111 const std::size_t N0 = diagonalizer->source()->matrixSize();
112 const std::size_t N = diagonalizer->matrixSize();
113 const std::size_t NN = N*N;
114
115 cmatrix work(N, N0, wrk); // matrix object for the workspace
116
117 cdiagonal gamma;
118
119 // in the beginning the P matrix is zero
120 std::fill_n(P.data(), NN, dcomplex(0.0));
121
122 std::exception_ptr error;
123
124 PLASK_OMP_PARALLEL_FOR_1
125 for (int l = 0; l < int(diagonalizer->lcount); ++l) {
126 try {
127 if (!error) diagonalizer->diagonalizeLayer(l);
128 } catch(...) {
129 error = std::current_exception();
130 }
131 }
132 if (error) std::rethrow_exception(error);
133
134 #ifdef OPENMP_FOUND
135 write_debug("{}: Entering into single region of reflection search", solver->getId());
136 #endif
137
138 // If we do not use emitting, we have to set field at the edge to 0 and the apply PML
139 if (!emitting) {
140 gamma = diagonalizer->Gamma(solver->stack[start]);
141 // Apply PML
142 // F(0) + B(0) = 0 ==> P(0) = -I
143 for (std::size_t i = 0; i < N; i++) {
144 dcomplex g = gamma[i] * solver->vpml.factor;
145 P(i,i) = - exp(-2. * I * g * solver->vpml.size); // P = phas * (-I) * phas
146 }
147 assert(!P.isnan());
148
149 // Shift matrix by `pmldist`
150 for (std::size_t i = 0; i < N; i++) phas[i] = exp(-I*gamma[i]*solver->vpml.dist);
151 assert(!phas.isnan());
152 mult_diagonal_by_matrix(phas, P); mult_matrix_by_diagonal(P, phas); // P = phas * P * phas
153 }
154
155 if (storeP == STORE_ALL) saveP(start);
156
157 for (std::size_t n = start; n != end; n += inc) {
158 gamma = diagonalizer->Gamma(solver->stack[n]);
159 assert(!gamma.isnan());
160
161 assert(!P.isnan());
162
163 if (n != start) {
164 double H = solver->vbounds->at(n) - solver->vbounds->at(n-1);
165 for (std::size_t i = 0; i < N; i++) phas[i] = exp(-I*gamma[i]*H);
166 assert(!phas.isnan());
167 mult_diagonal_by_matrix(phas, P); mult_matrix_by_diagonal(P, phas); // P = phas * P * phas
168 }
169
170 // Further calculations must be done only if the adjacent layers are not the same
171 if (solver->stack[n] != solver->stack[n+inc] || (emitting && n == start)) {
172 // temp = invTE(n+1)*TE(n) * [ phas*P*phas + I ]
173 assert(!diagonalizer->TE(solver->stack[n]).isnan());
174 assert(!diagonalizer->invTE(solver->stack[n]).isnan());
175 for (std::size_t i = 0, ii = 0; i < N; i++, ii += (N+1)) P[ii] += 1.; // P = P + I
176 if (solver->stack[n] != solver->stack[n+inc]) {
177 mult_matrix_by_matrix(diagonalizer->TE(solver->stack[n]), P, work); // work = TE[n] * P
178 mult_matrix_by_matrix(diagonalizer->invTE(solver->stack[n+inc]), work, temp);// temp = invTE[n+1] * work (= A)
179 } else {
180 std::copy_n(P.data(), NN, temp.data());
181 }
182
183 // P = invTH(n+1)*TH(n) * [ phas*P*phas - I ]
184 assert(!diagonalizer->TH(solver->stack[n]).isnan());
185 assert(!diagonalizer->invTH(solver->stack[n+inc]).isnan());
186 for (std::size_t i = 0, ii = 0; i < N; i++, ii += (N+1)) P[ii] -= 2.; // P = P - I
187
188 // multiply rows of P by -1 where necessary for properly outgoing wave
189 if (emitting && n == start) {
190 for (std::size_t i = 0; i < N; i++)
191 if (real(gamma[i]) < -SMALL)
192 for(std::size_t j = 0; j < N; j++) P(i,j) = -P(i,j);
193 }
194
195 if (solver->stack[n] != solver->stack[n+inc]) {
196 mult_matrix_by_matrix(diagonalizer->TH(solver->stack[n]), P, work); // work = TH[n] * P
197 mult_matrix_by_matrix(diagonalizer->invTH(solver->stack[n+inc]), work, P);// P = invTH[n+1] * work (= P)
198 }
199
200 // temp := temp-P, P := temp+P
201 for (std::size_t i = 0; i < NN; i++) {
202 dcomplex e = temp[i], h = P[i];
203 temp[i] = e - h;
204 P[i] = e + h;
205 }
206
207 // P = P * inv(temp)
208 int info;
209 zgetrf(int(N), int(N), temp.data(), int(N), ipiv, info); // temp = LU(temp)
210 if (info > 0) throw ComputationError(solver->getId(), "findReflection: Matrix [e(n) - h(n)] is singular");
211 assert(info == 0);
212 ztrsm('R', 'U', 'N', 'N', int(N), int(N), 1., temp.data(), int(N), P.data(), int(N)); // P = P * U^{-1}
213 ztrsm('R', 'L', 'N', 'U', int(N), int(N), 1., temp.data(), int(N), P.data(), int(N)); // P = P * L^{-1}
214 if (P.isnan()) throw ComputationError(solver->getId(), "findReflection: NaN in reflection matrix");
215 // reorder columns (there is no such function in LAPACK)
216 for (std::ptrdiff_t j = N-1; j >= 0; j--) {
217 int jp = ipiv[j]-1;
218 for (std::size_t i = 0; i < N; i++) std::swap(P(i,j), P(i,jp));
219 }
220 }
221
222 if (storeP == STORE_ALL) saveP(n+inc);
223 }
224 if (storeP == STORE_LAST) saveP(store);
225}
226
227
228 cvector ReflectionTransfer::getReflectionVector(const cvector& incident, IncidentDirection side)
229{
230 std::size_t last, first;
231
232 initDiagonalization();
233 switch (side) {
234 case INCIDENCE_TOP:
235 last = 0; first = solver->stack.size()-1; break;
236 case INCIDENCE_BOTTOM:
237 last = solver->stack.size()-1; first = 0; break;
238 }
239 findReflection(last, first, true);
240 return P * incident;
241}
242
243
244 cvector ReflectionTransfer::getTransmissionVector(const cvector& incident, IncidentDirection side)
245{
246 determineReflectedFields(incident, side);
247 ptrdiff_t n = (side == INCIDENCE_BOTTOM)? solver->stack.size()-1 : 0;
248 return
249 (((side == INCIDENCE_BOTTOM && n < solver->interface) ||
250 (side == INCIDENCE_TOP && n >= solver->interface))?
251 fields[n].F : fields[n].B);
252}
253
254
255// Some aliases
256#define F1 fields[n].F
257#define B1 fields[n].B
258#define F2 fields[n+inc].F
259#define B2 fields[n+inc].B
260
261
262 void ReflectionTransfer::determineFields()
263{
264 if (fields_determined == DETERMINED_RESONANT) return;
265
266 writelog(LOG_DETAIL, solver->getId() + ": Determining optical fields");
267
268 const std::size_t N = diagonalizer->matrixSize();
269 const std::size_t N0 = diagonalizer->source()->matrixSize();
270 const std::size_t NN = N*N;
271 cvector temp(wrk, N);
272
273 cdiagonal gamma;
274
275 size_t count = solver->stack.size();
276
277 // Assign the space for the field vectors
278 fields.resize(count);
279
280 // Obtain the physical fields at the last layer
281 storeP = STORE_LAST;
282 memP.resize(2);
283 interface_field = nullptr;
284 auto EH = getInterfaceVector();
285
286 for (unsigned pass = 0; pass < 1 || (pass < 2 && solver->interface != std::ptrdiff_t(count)); pass++)
287 {
288 // each pass for below and above the interface
289
290 std::size_t start, end;
291 std::ptrdiff_t inc;
292 switch (pass)
293 {
294 case 0: start = solver->interface-1; end = 0; inc = -1; break;
295 case 1: start = solver->interface; end = count-1; inc = +1; break;
296 }
297
298 fields[start].F.reset(N);
299 fields[start].B.reset(N);
300
301 // compute B-field for the layer next to the interface
302 std::size_t curr = solver->stack[start];
303
304 gamma = diagonalizer->Gamma(curr);
305
306 double H = (start == 0 || start == count-1)? 0. : (solver->vbounds->at(start) - solver->vbounds->at(start-1));
307 for (std::size_t i = 0; i < N; i++)
308 phas[i] = exp(-I*gamma[i]*H);
309
310 double II = (matching == MATCH_ADMITTANCE)? 1. : -1.;
311
312 // P = phas*P*phas + I
313 assert(memP[pass].rows() == N && memP[pass].cols() == N);
314 memcpy(P.data(), memP[pass].data(), NN*sizeof(dcomplex));
315 mult_diagonal_by_matrix(phas, P); mult_matrix_by_diagonal(P, phas); // P := phas * P * phas
316 for (std::size_t i = 0, ii = 0; i < N; i++, ii += (N+1)) P[ii] += II; // P := P ± I
317
318 if (matching == MATCH_ADMITTANCE)
319 mult_matrix_by_vector(diagonalizer->invTE(curr), EH, fields[start].B); // B := invTE * E
320 else
321 mult_matrix_by_vector(diagonalizer->invTH(curr), EH, fields[start].B); // B := invTH * H
322 invmult(P, fields[start].B); // B := inv(P) * B
323 for (std::size_t i = 0; i < N; i++) fields[start].B[i] *= phas[i]; // B := phas * B
324
325 // F-field for the first layer
326 mult_matrix_by_vector(memP[pass], fields[start].B, fields[start].F);
327
328 for (std::size_t n = start; n != end; n += inc) {
329 // Compute F and B field for the next (previous) layer
330
331 F2.reset(N);
332 B2.reset(N);
333
334 curr = solver->stack[n];
335 std::size_t next = solver->stack[n+inc];
336 assert(diagonalizer->isDiagonalized(curr));
337 assert(diagonalizer->isDiagonalized(next));
338
339 gamma = diagonalizer->Gamma(next);
340
341 if (next != curr) {
342 for (std::size_t i = 0; i < N; i++) F2[i] = F1[i] - B1[i]; // F2 := F1 - B1
343 mult_matrix_by_vector(diagonalizer->TH(curr), F2, temp); // temp := TH * F2
344 mult_matrix_by_vector(diagonalizer->invTH(next), temp, B2); // B2 := invTH * temp
345
346 for (std::size_t i = 0; i < N; i++) F2[i] = F1[i] + B1[i]; // F2 := F1 + B1
347 mult_matrix_by_vector(diagonalizer->TE(curr), F2, temp); // temp := TE * F2
348 memcpy(F2.data(), B2.data(), N*sizeof(dcomplex));
349 zgemm('N','N', int(N), 1, int(N0), 1., diagonalizer->invTE(next).data(), int(N),
350 temp.data(), int(N0), -1., B2.data(), int(N)); // B2 := invTE * temp - B2
351 for (std::size_t i = 0; i < N; i++)
352 F2[i] += 0.5 * B2[i]; // F2 := B2 + tH (F1-B2)
353 } else {
354 for (std::size_t i = 0; i < N; i++) B2[i] = 2. * B1[i];
355 memcpy(F2.data(), F1.data(), N*sizeof(dcomplex));
356 }
357
358 H = (n+inc == end)? 0. : (solver->vbounds->at(n+inc) - solver->vbounds->at(n+inc-1));
359 for (std::size_t i = 0; i < N; i++) {
360 dcomplex phas = exp(-I*gamma[i]*H);
361 B2[i] *= 0.5 * phas; // B2 := 1/2 * phas * B2
362 F2[i] /= phas; // F2 := phas^(-1) * F2
363 }
364 }
365 }
366
367 storeP = STORE_NONE;
368 fields_determined = DETERMINED_RESONANT;
369
370 // Finally normalize fields
371 if (solver->emission == ModalBase::EMISSION_BOTTOM || solver->emission == ModalBase::EMISSION_TOP) {
372 size_t n = (solver->emission == ModalBase::EMISSION_BOTTOM)? 0 : count-1;
373
374 double P = 1./Z0 * abs(diagonalizer->source()->integratePoyntingVert(getFieldVectorE(0., n),
375 getFieldVectorH(0., n)));
376
377 if (P < SMALL) {
378 writelog(LOG_WARNING, "Device is not emitting to the {} side: skipping normalization",
379 (solver->emission == ModalBase::EMISSION_TOP)? "top" : "bottom");
380 } else {
381 P = 1. / sqrt(P);
382 for (size_t n = 0; n < count; ++n) {
383 F1 *= P;
384 B1 *= P;
385 }
386 }
387 }
388}
389
390
391 void ReflectionTransfer::determineReflectedFields(const cvector& incident, IncidentDirection side)
392{
393 if (fields_determined == DETERMINED_REFLECTED && incident == incident_vector && side == incident_side) return;
394 incident_vector = incident.copy();
395 incident_side = side;
396
397 writelog(LOG_DETAIL, solver->getId() + ": Determining reflected optical fields");
398
399 size_t count = solver->stack.size();
400
401 // Assign the space for the field vectors
402 fields.resize(count);
403
404 std::size_t start, end;
405 std::ptrdiff_t inc;
406 switch (side)
407 {
408 case INCIDENCE_TOP: start = count-1; end = 0; inc = -1; break;
409 case INCIDENCE_BOTTOM: start = 0; end = count-1; inc = +1; break;
410 }
411
412 // Store all reflectivities
413 storeP = STORE_ALL;
414 memP.resize(count);
415
416 // Compute reflection matrices
417 initDiagonalization();
418 findReflection(end, start, true);
419
420 // Temporary and initial data
421 const std::size_t N = diagonalizer->matrixSize();
422 const std::size_t N0 = diagonalizer->source()->matrixSize();
423 cvector temp(wrk, N);
424 cdiagonal gamma;
425
426 std::size_t curr = solver->stack[start];
427 double H;
428
429 fields[start].B = incident.copy(); // diagonalized incident E-field
430 fields[start].F.reset(N);
431
432 for (std::size_t n = start; n != end; n += inc)
433 {
434 // F-field for the current layer
435 mult_matrix_by_vector(memP[n], B1, F1);
436
437 // Compute B field for the next (previous) layer
438
439 F2.reset(N);
440 B2.reset(N);
441
442 curr = solver->stack[n];
443 const std::size_t next = solver->stack[n+inc];
444
445 gamma = diagonalizer->Gamma(next);
446
447 if (next != curr || n+inc == end) {
448 if (next != curr) {
449 for (std::size_t i = 0; i < N; i++) F2[i] = F1[i] - B1[i]; // F2 := F1 - B1
450 mult_matrix_by_vector(diagonalizer->TH(curr), F2, temp); // temp := TH * F2
451 mult_matrix_by_vector(diagonalizer->invTH(next), temp, B2); // B2 := invTH * temp
452 } else {
453 for (std::size_t i = 0; i < N; i++) B2[i] = F1[i] - B1[i]; // B2 := F1 - B1
454 }
455 // multiply rows of invTH by -1 where necessary for the outer layer
456 if (n+inc == end) {
457 for (std::size_t i = 0; i < N; i++)
458 if (real(gamma[i]) < -SMALL) B2[i] = -B2[i];
459 }
460
461 for (std::size_t i = 0; i < N; i++) F2[i] = F1[i] + B1[i]; // F2 := F1 + B1
462 if (next != curr) {
463 mult_matrix_by_vector(diagonalizer->TE(curr), F2, temp); // temp := TE * F2
464 zgemm('N','N', int(N), 1, int(N0), 1., diagonalizer->invTE(next).data(), int(N),
465 temp.data(), int(N0), -1., B2.data(), int(N)); // B2 := invTE * temp - B2
466 } else {
467 for (std::size_t i = 0; i < N; i++)
468 B2[i] = F2[i] - B2[i]; // B2 := (F1+B1) + (F1-B1)
469 }
470 } else {
471 for (std::size_t i = 0; i < N; i++) B2[i] = 2. * B1[i];
472 }
473
474 if (n+inc != end) {
475 H = solver->vbounds->at(n+inc) - solver->vbounds->at(n+inc-1);
476 for (std::size_t i = 0; i < N; i++)
477 B2[i] *= 0.5 * exp(-I*gamma[i]*H); // B2 := 1/2 * phas * B2
478 } else {
479 for (std::size_t i = 0; i < N; i++) B2[i] *= 0.5; // B2 := 1/2 * phas * B2
480 }
481 }
482
483 fields[end].F = cvector(N, 0.);
484
485 // In the outer layers replace F and B where necessary for consistent gamma handling
486 for (std::size_t n = 0; n < count; n += count-1) {
487 gamma = diagonalizer->Gamma(solver->stack[n]);
488 for (std::size_t i = 0; i < N; i++) {
489 if (real(gamma[i]) < -SMALL)
490 std::swap(fields[n].F, fields[n].B);
491 }
492 }
493
494 // Replace F and B at one side of the interface for consistency in getFieldVectorE and getFieldVectorH
495 size_t interface = size_t(max(solver->interface, ptrdiff_t(0)));
496 switch (side)
497 {
498 case INCIDENCE_TOP: start = interface; end = count; break;
499 case INCIDENCE_BOTTOM: start = 0; end = min(interface, count); break;
500 }
501 // start = size_t(max(solver->interface, ptrdiff_t(0))); end = count;
502 for (std::size_t n = start; n < end; n++) {
503 gamma = diagonalizer->Gamma(solver->stack[n]);
504 H = (n < count-1 && n > 0)? solver->vbounds->at(n) - solver->vbounds->at(n-1) : 0.;
505 for (std::size_t i = 0; i < N; i++) {
506 dcomplex phas = exp(-I*gamma[i]*H);
507 dcomplex t = B1[i] / phas;
508 if (isnan(t) && B1[i] == 0.) t = 0.;
509 B1[i] = F1[i] * phas;
510 if (isnan(B1[i]) && F1[i] == 0.) B1[i] = 0.;
511 F1[i] = t;
512 }
513 }
514
515 storeP = STORE_NONE;
516 fields_determined = DETERMINED_REFLECTED;
517}
518
519
520 cvector ReflectionTransfer::getFieldVectorE(double z, std::size_t n, PropagationDirection part)
521{
522 assert(fields_determined != DETERMINED_NOTHING);
523
524 adjust_z(n, z, part);
525
526 cdiagonal gamma = diagonalizer->Gamma(solver->stack[n]);
527
528 const std::size_t N = gamma.size();
529 cvector E(N);
530
531 for (std::size_t i = 0; i < N; i++) {
532 dcomplex phi = - I * gamma[i] * z;
533 dcomplex ef = F1[i] * exp(phi), eb = B1[i] * exp(-phi);
534 if (part == PROPAGATION_UPWARDS) eb = 0;
535 else if (part == PROPAGATION_DOWNWARDS) ef = 0;
536 if (isnan(ef) || isinf(ef.real()) || isinf(ef.imag())) ef = exp(log(F1[i]) + phi);
537 if (isnan(eb) || isinf(eb.real()) || isinf(eb.imag())) eb = exp(log(B1[i]) - phi);
538 if (isnan(ef) || isinf(ef.real()) || isinf(ef.imag())) ef = 0.; // not elegant but allows to avoid NaNs
539 if (isnan(eb) || isinf(eb.real()) || isinf(eb.imag())) eb = 0.; // not elegant but allows to avoid NaNs
540 E[i] = ef + eb;
541 }
542
543 return diagonalizer->TE(solver->stack[n]) * E;
544}
545
546
547 cvector ReflectionTransfer::getFieldVectorH(double z, std::size_t n, PropagationDirection part)
548{
549 assert(fields_determined != DETERMINED_NOTHING);
550
551 adjust_z(n, z, part);
552
553 cdiagonal gamma = diagonalizer->Gamma(solver->stack[n]);
554
555 const std::size_t N = gamma.size();
556 cvector H(N);
557
558 for (std::size_t i = 0; i < N; i++) {
559 dcomplex phi = - I * gamma[i] * z;
560 dcomplex ef = F1[i] * exp(phi), eb = B1[i] * exp(-phi);
561 if (part == PROPAGATION_UPWARDS) eb = 0;
562 else if (part == PROPAGATION_DOWNWARDS) ef = 0;
563 if (isnan(ef) || isinf(ef.real()) || isinf(ef.imag())) ef = exp(log(F1[i]) + phi);
564 if (isnan(eb) || isinf(eb.real()) || isinf(eb.imag())) eb = exp(log(B1[i]) - phi);
565 if (isnan(ef) || isinf(ef.real()) || isinf(ef.imag())) ef = 0.; // not elegant but allows to avoid NaNs
566 if (isnan(eb) || isinf(eb.real()) || isinf(eb.imag())) eb = 0.; // not elegant but allows to avoid NaNs
567 H[i] = ef - eb;
568 }
569
570 if (n == 0 || std::size_t(n) == solver->vbounds->size()) {
571 // In the outer layers multiply H by -1 where necessary for propagating wave
572 for (std::size_t i = 0; i < N; i++)
573 if (real(gamma[i]) < -SMALL) H[i] = - H[i];
574 }
575
576 return diagonalizer->TH(solver->stack[n]) * H;
577}
578
579
580 inline static dcomplex _int_refl(double z1, double z2, dcomplex g, dcomplex E) {
581 if (is_zero(E)) return 0.;
582 if (is_zero(g)) return E * (z2 - z1);
583 dcomplex res = -I * E / g * (exp(I * g * z2) - exp(I * g * z1));
584 if (isinf(res.real()) || isinf(res.imag()) || isnan(res)) {
585 dcomplex logE = log(E);
586 res = -I / g * (exp(logE + I * g * z2) - exp(logE + I * g * z1));
587 }
588 return res;
589}
590
591 double ReflectionTransfer::integrateField(WhichField field, size_t n, double z1, double z2) {
592 size_t layer = solver->stack[n];
593 size_t N = diagonalizer->matrixSize();
594
595 cmatrix TE = diagonalizer->TE(layer),
596 TH = diagonalizer->TH(layer);
597 cdiagonal gamma = diagonalizer->Gamma(layer);
598 assert(gamma.size() == N);
599
600 adjust_z(n, z1, z2);
601
602 return diagonalizer->source()->integrateField(field, layer, TE, TH,
603 [n, z1, z2, gamma, this](size_t i, size_t j) {
604 return std::make_pair(
605 _int_refl(z1, z2, -gamma[i] + conj(gamma[j]), F1[i] * conj(F1[j])) +
606 _int_refl(z1, z2, gamma[i] + conj(gamma[j]), B1[i] * conj(F1[j])) +
607 _int_refl(z1, z2, -gamma[i] - conj(gamma[j]), F1[i] * conj(B1[j])) +
608 _int_refl(z1, z2, gamma[i] - conj(gamma[j]), B1[i] * conj(B1[j])),
609
610 _int_refl(z1, z2, -gamma[i] + conj(gamma[j]), F1[i] * conj( F1[j])) +
611 _int_refl(z1, z2, gamma[i] + conj(gamma[j]), -B1[i] * conj( F1[j])) +
612 _int_refl(z1, z2, -gamma[i] - conj(gamma[j]), F1[i] * conj(-B1[j])) +
613 _int_refl(z1, z2, gamma[i] - conj(gamma[j]), -B1[i] * conj(-B1[j]))
614 );
615 });
616}
617
618}}} // namespace plask::optical::modal