PLaSK library

Loading...
Searching...
No Matches
toeplitz.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_TOEPLITZ_H
15#define PLASK__SOLVER_OPTICAL_MODAL_TOEPLITZ_H
16
17#include <plask/plask.hpp>
18#include "../matrices.hpp"
19
20 namespace plask { namespace optical { namespace modal {
21
28 template <typename D>
29 void ToeplitzLevinson(const DataVector<D>& TM, Matrix<D>& X)
30{
31 assert(TM.size() % 2 == 1);
32 const size_t n = (TM.size() + 1) / 2;
33 const size_t n2 = TM.size() - 1;
34 const size_t nx = X.cols();
35
36 std::unique_ptr<D[]> F(new D[n]);
37 std::unique_ptr<D[]> B(new D[n]);
38 std::unique_ptr<D[]> eX(new D[nx]);
39
40 if (TM[0] == 0.) throw ComputationError("ToeplitzLevinson", "cannot invert Fourier coefficients matrix");
41
42 F[0] = 1. / TM[0];
43 B[0] = 1. / TM[0];
44 for (size_t r = 0, dx = 0; r < nx; ++r, dx += n)
45 X[dx] /= TM[0];
46
47 for (size_t i = 1; i < n; i++) {
48 F[i] = B[i] = 0.;
49 D ef = 0., eb = 0.;
50
51 for (size_t r = 0; r < nx; r++) eX[r] = 0.;
52
53 for (size_t j = 0; j < i; j++) {
54 size_t ij = i-j;
55 ef += TM[ij] * F[j];
56 eb += TM[n2-j] * B[j];
57 for (size_t r = 0, mj = j; r < nx; ++r, mj += n)
58 eX[r] += TM[ij] * X[mj];
59 }
60
61 D scal = (1. - ef * eb);
62 if (scal == 0.) throw ComputationError("ToeplitzLevinson", "cannot invert Fourier coefficients matrix");
63 scal = 1. / scal;
64
65 D b = B[0];
66 B[0] = -(eb * scal) * F[0];
67 F[0] *= scal;
68 for (size_t j = 1; j <= i; j++) {
69 D f = F[j];
70 D bj = B[j];
71 F[j] = (scal * f) - ((ef * scal) * b);
72 B[j] = (scal * b) - ((eb * scal) * f);
73 b = bj;
74 }
75
76 for (size_t r = 0, dx = 0; r < nx; ++r, dx += n) {
77 size_t ix = dx + i;
78 for (size_t j = 0; j < i; j++)
79 X[dx+j] += (X[ix] - eX[r]) * B[j];
80 X[ix] = (X[ix] - eX[r]) * B[i];
81 }
82 }
83}
84
85
86
87}}} // namespace plask::optical::modal
88#endif // PLASK__SOLVER_OPTICAL_MODAL_TOEPLITZ_H