PLaSK library
Loading...
Searching...
No Matches
expansioncyl-infini.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#include "expansioncyl-infini.hpp"
15#include "solvercyl.hpp"
16#include "zeros-data.hpp"
17#include "../plask/math.hpp"
18
19#include "../gauss_legendre.hpp"
20#include "../gauss_laguerre.hpp"
21
22#include "besselj.hpp"
23
24#define SOLVER static_cast<BesselSolverCyl*>(solver)
25
26namespace plask { namespace optical { namespace modal {
27
31
32
33void ExpansionBesselInfini::init2()
34{
35 SOLVER->writelog(LOG_DETAIL, "Preparing Bessel functions for m = {}", m);
36
37 if (SOLVER->geometry->getEdge(Geometry::DIRECTION_TRAN, true).type() != edge::Strategy::DEFAULT &&
38 SOLVER->geometry->getEdge(Geometry::DIRECTION_TRAN, true).type() != edge::Strategy::SIMPLE &&
39 SOLVER->geometry->getEdge(Geometry::DIRECTION_TRAN, true).type() != edge::Strategy::EXTEND)
40 throw BadInput(solver->getId(), "outer geometry edge must be 'extend' or a simple material");
41
42 double k0 = isnan(lam0)? this->k0.real() : 2e3*M_PI / lam0;
43 double kmax = SOLVER->kmax * k0;
44
45 size_t N = SOLVER->size;
46 double R = rbounds[rbounds.size()-1];
47 double ib = 1. / R;
48 double kdlt;
49
50 switch (SOLVER->kmethod) {
51 case BesselSolverCyl::WAVEVECTORS_UNIFORM:
52 SOLVER->writelog(LOG_DETAIL, "Using uniform wavevectors");
53 if (isnan(k0)) throw BadInput(SOLVER->getId(), "no wavelength given: specify 'lam' or 'lam0'");
54 kpts.resize(N);
55 kdlt = SOLVER->kscale * kmax / N;
56 kdelts.reset(N, kdlt);
57 kdlt *= R;
58 for (size_t i = 0; i != N; ++i) kpts[i] = (0.5 + double(i) * kdlt);
59 break;
60 case BesselSolverCyl::WAVEVECTORS_LAGUERRE:
61 SOLVER->writelog(LOG_DETAIL, "Using Laguerre wavevectors");
62 gaussLaguerre(N, kpts, kdelts, 1. / (SOLVER->kscale));
63 kdelts *= ib;
64 break;
65 case BesselSolverCyl::WAVEVECTORS_MANUAL:
66 SOLVER->writelog(LOG_DETAIL, "Using manual wavevectors");
67 kpts.resize(N);
68 kdelts.reset(N);
69 if (!SOLVER->kweights) {
70 if (SOLVER->klist.size() != N+1)
71 throw BadInput(SOLVER->getId(), "if no weights are given, number of manually specified wavevectors must be {}",
72 N+1);
73 for (size_t i = 0; i != N; ++i) {
74 kpts[i] = 0.5 * (SOLVER->klist[i] + SOLVER->klist[i+1]) * R;
75 kdelts[i] = (SOLVER->klist[i+1] - SOLVER->klist[i]);
76 }
77 } else {
78 if (SOLVER->klist.size() != N)
79 throw BadInput(SOLVER->getId(), "if weights are given, number of manually specified wavevectors must be {}",
80 N);
81 if (SOLVER->kweights->size() != N)
82 throw BadInput(SOLVER->getId(), "number of manually specified wavevector weights must be {}", N+1);
83 kpts = SOLVER->klist;
84 for (double& k: kpts) k *= R;
85 kdelts.reset(SOLVER->kweights->begin(), SOLVER->kweights->end());
86 }
87 break;
88 case BesselSolverCyl::WAVEVECTORS_NONUNIFORM:
89 SOLVER->writelog(LOG_DETAIL, "Using non-uniform wavevectors");
90 if (isnan(k0)) throw BadInput(SOLVER->getId(), "no wavelength given: specify 'lam' or 'lam0'");
91 kpts.resize(N);
92 kdelts.reset(N);
93 // Häyrynen, T., de Lasson, J.R., Gregersen, N., 2016.
94 // Open-geometry Fourier modal method: modeling nanophotonic structures in infinite domains.
95 // J. Opt. Soc. Am. A 33, 1298. https://doi.org/10.1364/josaa.33.001298
96 int M1, M2, M3;
97 M1 = M2 = (N+1) / 3;
98 M3 = N - M1 - M2;
99 if (M3 < 0) throw BadInput(SOLVER->getId(), "too small expansion size");
100 int i = 0;
101 double k1 = 0.;
102 for(int m = 1; m <= M1; ++m, ++i) {
103 double theta = 0.5*M_PI * double(m) / M1;
104 double k2 = k0 * sin(theta) * SOLVER->kscale;
105 kdelts[i] = k2 - k1;
106 kpts[i] = 0.5 * (k1 + k2) * R;
107 k1 = k2;
108 }
109 double dlt1;
110 for(int m = 1; m <= M2; ++m, ++i) {
111 double theta = 0.5*M_PI * (1. + double(m) / M2);
112 double k2 = k0 * (2 - sin(theta)) * SOLVER->kscale;
113 kdelts[i] = dlt1 = k2 - k1;
114 kpts[i] = 0.5 * (k1 + k2) * R;
115 k1 = k2;
116 }
117 double km = k1;
118 double dlt2 = (kmax * SOLVER->kscale - km - M3 * dlt1) / (M3 * (M3 + 1));
119 if (dlt2 < 0)
120 throw BadInput(SOLVER->getId(), "for non-uniform wavevectors kmax must be at least {}",
121 (km + M3 * dlt1) / SOLVER->kscale / k0);
122 for(int m = 1; m <= M3; ++m, ++i) {
123 double k2 = km + dlt1 * m + dlt2 * m * (m+1);
124 kdelts[i] = k2 - k1;
125 kpts[i] = 0.5 * (k1 + k2) * R;
126 k1 = k2;
127 }
128 break;
129 }
130
131 init3();
132}
133
134
135void ExpansionBesselInfini::getMatrices(size_t layer, cmatrix& RE, cmatrix& RH)
136{
137 assert(initialized);
138 if (isnan(k0)) throw BadInput(SOLVER->getId(), "wavelength or k0 not set");
139 if (isinf(k0.real())) throw BadInput(SOLVER->getId(), "wavelength must not be 0");
140
141 size_t N = SOLVER->size;
142 dcomplex ik0 = 1. / k0;
143 double ib = 1. / rbounds[rbounds.size()-1];
144
145 const Integrals& eps = layers_integrals[layer];
146
147 for (size_t j = 0; j != N; ++j) {
148 size_t js = idxs(j), jp = idxp(j);
149 // double k = kpts[j] * ib;
150 for (size_t i = 0; i != N; ++i) {
151 size_t is = idxs(i), ip = idxp(i);
152 double g = kpts[i] * ib;
153 dcomplex gVk = 0.5 * ik0 * g * eps.V_k(i,j);
154 RH(is,jp) = gVk;
155 RH(is,js) = gVk;
156 RH(ip,jp) = -gVk;
157 RH(ip,js) = -gVk;
158 }
159 RH(js,js) -= k0;
160 RH(jp,jp) += k0;
161 }
162
163 for (size_t j = 0; j != N; ++j) {
164 size_t js = idxs(j), jp = idxp(j);
165 for (size_t i = 0; i != N; ++i) {
166 size_t is = idxs(i), ip = idxp(i);
167 // double g = kpts[i] * ib;
168 RE(ip,js) = 0.5 * k0 * eps.Tps(i,j);
169 RE(ip,jp) = 0.5 * k0 * eps.Tpp(i,j);
170 RE(is,js) = -0.5 * k0 * eps.Tss(i,j);
171 RE(is,jp) = -0.5 * k0 * eps.Tsp(i,j);
172 }
173 double k = kpts[j] * ib;
174 dcomplex ik0k2 = 0.5 * ik0 * k*k;
175 RE(jp,js) -= ik0k2;
176 RE(jp,jp) -= ik0k2;
177 RE(js,js) += ik0k2;
178 RE(js,jp) += ik0k2;
179 }
180}
181
182void ExpansionBesselInfini::integrateParams(Integrals& integrals,
183 const dcomplex* datap, const dcomplex* datar, const dcomplex* dataz,
184 dcomplex datap0, dcomplex datar0, dcomplex dataz0) {
185 auto raxis = mesh->tran();
186
187 size_t nr = raxis->size(), N = SOLVER->size;
188 double R = rbounds[rbounds.size()-1];
189 double ib = 1. / R;
190 double* factors; // scale factors for making matrices orthonormal
191
192 integrals.reset(N);
193
194 TempMatrix temp = getTempMatrix();
195 aligned_unique_ptr<double> _tmp;
196
197 cmatrix JmJp, JpJm;
198 if (N < 2) {
199 _tmp.reset(aligned_malloc<double>(4*N));
200 factors = _tmp.get();
201 } else if (SOLVER->rule == BesselSolverCyl::RULE_COMBINED_1) {
202 factors = reinterpret_cast<double*>(integrals.V_k.data());
203 } else if (SOLVER->rule == BesselSolverCyl::RULE_COMBINED_2) {
204 _tmp.reset(aligned_malloc<double>(3*N + 4*N*N));
205 JmJp.reset(N, N, integrals.V_k.data());
206 JpJm.reset(N, N, reinterpret_cast<dcomplex*>(_tmp.get() + 3*N));
207 factors = _tmp.get();
208 } else {
209 factors = reinterpret_cast<double*>(temp.data());
210 }
211 for (size_t i = 0; i < N; ++i) {
212 factors[i] = kpts[i] * ib * kdelts[i];
213 }
214 double* Jm = factors + N;
215 double* Jp = factors + 2*N;
216
217 if (SOLVER->rule == BesselSolverCyl::RULE_OLD) {
218
219 double* J = factors + 3*N;
220
221 zero_matrix(integrals.V_k);
222 zero_matrix(integrals.Tss);
223 zero_matrix(integrals.Tsp);
224 zero_matrix(integrals.Tps);
225 zero_matrix(integrals.Tpp);
226
227 for (size_t ri = 0; ri != nr; ++ri) {
228 double r = raxis->at(ri);
229 dcomplex repst = r * (datar[ri] + datap[ri]), repsd = r * (datar[ri] - datap[ri]);
230 const dcomplex riepsz = r * dataz[ri];
231
232 for (size_t i = 0; i < N; ++i) {
233 double kr = kpts[i] * ib * r;
234 Jm[i] = cyl_bessel_j(m-1, kr);
235 J[i] = cyl_bessel_j(m, kr);
236 Jp[i] = cyl_bessel_j(m+1, kr);
237 }
238 for (size_t j = 0; j < N; ++j) {
239 double k = kpts[j] * ib;
240 double Jk = J[j], Jmk = Jm[j], Jpk = Jp[j];
241 for (size_t i = 0; i < N; ++i) {
242 double Jg = J[i], Jmg = Jm[i], Jpg = Jp[i];
243 integrals.V_k(i,j) += Jg * riepsz * Jk * k * factors[i];
244 integrals.Tss(i,j) += Jmg * repst * Jmk * factors[i];
245 integrals.Tsp(i,j) += Jmg * repsd * Jpk * factors[i];
246 integrals.Tps(i,j) += Jpg * repsd * Jmk * factors[i];
247 integrals.Tpp(i,j) += Jpg * repst * Jpk * factors[i];
248 }
249 }
250 }
251 for (size_t i = 0; i < N; ++i) {
252 integrals.V_k(i,i) += dataz0 * kpts[i] * ib;
253 dcomplex epst = datar0 + datap0, epsd = datar0 - datap0;
254 integrals.Tss(i,i) += epst;
255 integrals.Tsp(i,i) += epsd;
256 integrals.Tps(i,i) += epsd;
257 integrals.Tpp(i,i) += epst;
258 }
259
260 } else {
261
262 if (SOLVER->rule == BesselSolverCyl::RULE_DIRECT) {
263
264 zero_matrix(integrals.Tss);
265 zero_matrix(integrals.Tsp);
266 zero_matrix(integrals.Tps);
267 zero_matrix(integrals.Tpp);
268
269 for (size_t ri = 0; ri != nr; ++ri) {
270 double r = raxis->at(ri);
271 dcomplex repst = r * (datar[ri] + datap[ri]), repsd = r * (datar[ri] - datap[ri]);
272 for (size_t i = 0; i < N; ++i) {
273 double kr = kpts[i] * ib * r;
274 Jm[i] = cyl_bessel_j(m-1, kr);
275 Jp[i] = cyl_bessel_j(m+1, kr);
276 }
277 for (size_t i = 0; i < N; ++i) {
278 dcomplex cs = factors[i] * repst, cd = factors[i] * repsd;
279 for (size_t j = 0; j < N; ++j) {
280 integrals.Tss(i,j) += cs * Jm[i] * Jm[j];
281 integrals.Tsp(i,j) += cd * Jm[i] * Jp[j];
282 integrals.Tps(i,j) += cd * Jp[i] * Jm[j];
283 integrals.Tpp(i,j) += cs * Jp[i] * Jp[j];
284 }
285 }
286 }
287 for (size_t i = 0; i < N; ++i) {
288 dcomplex epst = datar0 + datap0, epsd = datar0 - datap0;
289 integrals.Tss(i,i) += epst;
290 integrals.Tsp(i,i) += epsd;
291 integrals.Tps(i,i) += epsd;
292 integrals.Tpp(i,i) += epst;
293 }
294
295 } else {
296
297 if (SOLVER->rule == BesselSolverCyl::RULE_COMBINED_1) {
298
299 cmatrix workess(N, N, temp.data()), workepp(N, N, temp.data()+N*N),
300 worksp(N, N, temp.data()+2*N*N), workps(N, N, temp.data()+3*N*N);
301
302 zero_matrix(workess);
303 zero_matrix(workepp);
304 zero_matrix(worksp);
305 zero_matrix(workps);
306
307 for (size_t ri = 0; ri != nr; ++ri) {
308 double r = raxis->at(ri);
309 dcomplex riepsr = r * datar[ri];
310 for (size_t i = 0; i < N; ++i) {
311 double kr = kpts[i] * ib * r;
312 Jm[i] = cyl_bessel_j(m-1, kr);
313 Jp[i] = cyl_bessel_j(m+1, kr);
314 }
315 for (size_t i = 0; i < N; ++i) {
316 dcomplex ce = factors[i] * riepsr;
317 for (size_t j = 0; j < N; ++j) {
318 workess(i,j) += ce * Jm[i] * Jm[j];
319 workepp(i,j) += ce * Jp[i] * Jp[j];
320 }
321 }
322 }
323 for (size_t i = 0; i < N; ++i) {
324 workess(i,i) += datar0;
325 workepp(i,i) += datar0;
326 }
327
328 make_unit_matrix(integrals.Tss);
329 make_unit_matrix(integrals.Tpp);
330
331 invmult(workess, integrals.Tss);
332 invmult(workepp, integrals.Tpp);
333
334 zero_matrix(integrals.Tsp);
335 zero_matrix(integrals.Tps);
336
337 if (N > 2) {
338 std::copy_n(factors, N, reinterpret_cast<double*>(temp.data()));
339 factors = reinterpret_cast<double*>(temp.data());
340 Jm = factors + N; Jp = factors + 2*N;
341 }
342
343 } else { // if (SOLVER->rule == BesselSolverCyl::RULE_COMBINED_2)
344
345 cmatrix work(temp);
346
347 zero_matrix(integrals.TT);
348 zero_matrix(work);
349 zero_matrix(JmJp);
350 zero_matrix(JpJm);
351
352 for (size_t ri = 0; ri != nr; ++ri) {
353 double r = raxis->at(ri);
354 dcomplex riepsr = r * datar[ri];
355 for (size_t i = 0; i < N; ++i) {
356 double kr = kpts[i] * ib * r;
357 Jm[i] = cyl_bessel_j(m-1, kr);
358 Jp[i] = cyl_bessel_j(m+1, kr);
359 }
360 for (size_t j = 0; j < N; ++j) {
361 for (size_t i = 0; i < N; ++i) {
362 dcomplex ce = riepsr * factors[i];
363 integrals.TT(i,j) += ce * Jm[i] * Jm[j];
364 integrals.TT(i,j+N) += ce * Jm[i] * Jp[j];
365 integrals.TT(i+N,j) += ce * Jp[i] * Jm[j];
366 integrals.TT(i+N,j+N) += ce * Jp[i] * Jp[j];
367 }
368 }
369 }
370 for (size_t i = 0; i < 2*N; ++i) {
371 integrals.TT(i,i) += datar0;
372 }
373
374 // Compute Jp(kr) Jm(gr) r dr and Jp(gr) Jm(kr) r dr using analytical formula
375 for (size_t j = 0; j < N; ++j) {
376 double k = kpts[j] * ib;
377 for (size_t i = 0; i < j; ++i) {
378 double g = kpts[i] * ib;
379 double val = factors[i] * 2*m / (k*g) * pow(g/k, m);
380 JmJp(i,j) = work(i,j+N) = val; // g<k s=g p=k
381 integrals.TT(i,j+N) += datar0 * val;
382 }
383 double val = factors[j] * m / (k*k) - 1.;
384 JmJp(j,j) = JpJm(j,j) = work(j,j+N) = work(j+N,j) = val;
385 dcomplex iepsr0 = val * datar0;
386 integrals.TT(j,j+N) += iepsr0;
387 integrals.TT(j+N,j) += iepsr0;
388 for (size_t i = j+1; i < N; ++i) {
389 double g = kpts[i] * ib;
390 double val = factors[i] * 2*m / (k*g) * pow(k/g, m);
391 JpJm(i,j) = work(i+N,j) = val; // k<g s=k p=g
392 integrals.TT(i+N,j) += datar0 * val;
393 }
394 }
395 for (size_t i = 0; i < N; ++i) work(i,i) = 1.;
396 for (size_t i = 0; i < N; ++i) work(i+N,i+N) = 1.;
397
398 invmult(integrals.TT, work);
399
400 zero_matrix(integrals.Tsp);
401 zero_matrix(integrals.Tps);
402
403 for (size_t j = 0; j < N; ++j) {
404 for (size_t i = 0; i < N; ++i) integrals.Tss(i,j) = work(i,j);
405 for (size_t i = 0; i < N; ++i) integrals.Tps(i,j) = work(i+N,j);
406 }
407
408 zgemm('N', 'N', int(N), int(N), int(N), 1., JpJm.data(), int(N), work.data()+2*N*N, int(2*N), 1.,
409 integrals.Tpp.data(), int(N));
410 zgemm('N', 'N', int(N), int(N), int(N), 1., JmJp.data(), int(N), work.data()+N, int(2*N), 1.,
411 integrals.Tss.data(), int(N));
412 }
413
414 for (size_t ri = 0; ri != nr; ++ri) {
415 double r = raxis->at(ri);
416 dcomplex repsp = r * datap[ri];
417 for (size_t i = 0; i < N; ++i) {
418 double kr = kpts[i] * ib * r;
419 Jm[i] = cyl_bessel_j(m-1, kr);
420 Jp[i] = cyl_bessel_j(m+1, kr);
421 }
422 for (size_t i = 0; i < N; ++i) {
423 dcomplex c = repsp * factors[i];
424 for (size_t j = 0; j < N; ++j) {
425 integrals.Tss(i,j) += c * Jm[i] * Jm[j];
426 integrals.Tpp(i,j) += c * Jp[i] * Jp[j];
427 }
428 }
429 }
430 for (size_t i = 0; i < N; ++i) {
431 integrals.Tss(i,i) += datap0;
432 integrals.Tpp(i,i) += datap0;
433 }
434 }
435
436 cmatrix work(N, N, temp.data()+N*N);
437
438 zero_matrix(work);
439 double* J = Jm;
440
441 for (size_t ri = 0; ri != nr; ++ri) {
442 double r = raxis->at(ri);
443 dcomplex repsz = r * dataz[ri];
444 for (size_t i = 0; i < N; ++i) {
445 double kr = kpts[i] * ib * r;
446 J[i] = cyl_bessel_j(m, kr);
447 }
448 for (size_t j = 0; j < N; ++j) {
449 for (size_t i = 0; i < N; ++i) {
450 work(i,j) += repsz * factors[i] * J[i] * J[j];
451 }
452 }
453 }
454 for (size_t i = 0; i < N; ++i) work(i,i) += dataz0;
455
456 // make_unit_matrix(integrals.V_k);
457 zero_matrix(integrals.V_k);
458 for (int i = 0; i < N; ++i) {
459 double k = kpts[i] * ib;
460 integrals.V_k(i,i) = k;
461 }
462 invmult(work, integrals.V_k);
463 // for (size_t i = 0; i < N; ++i) {
464 // double g = kpts[i] * ib;
465 // for (size_t j = 0; j < N; ++j) integrals.V_k(i,j) *= g;
466 // }
467
468 }
469}
470
471}}} // # namespace plask::optical::modal