PLaSK library

Loading...
Searching...
No Matches
gauss_laguerre.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 "gauss_laguerre.hpp"
15#include "fortran.hpp"
16
17#include <memory>
18
19 namespace plask { namespace optical { namespace modal {
20
21 inline static double sign(double x) { return (x < 0.)? -1. : 1.; }
22
79 inline static void imtqlx(int n, double d[], double e[], double z[])
80{
81 double b;
82 double c;
83 double f;
84 double g;
85 int i;
86 int ii;
87 const int itn = 30;
88 int j;
89 int k;
90 int l;
91 int m;
92 int mml;
93 double p;
94 double prec;
95 double r;
96 double s;
97
98 prec = std::numeric_limits<double>::epsilon();
99
100 if (n == 1) return;
101
102 e[n-1] = 0.0;
103
104 for (l = 1; l <= n; l++) {
105 j = 0;
106 for (;;) {
107 for (m = l; m <= n; m++) {
108 if (m == n) break;
109 if (fabs(e[m-1]) <= prec * (fabs(d[m-1]) + fabs(d[m])))
110 break;
111 }
112 p = d[l-1];
113 if (m == l) break;
114 if (j >= itn)
115 throw ComputationError("imtqlx", "iteration limit exceeded\n");
116 j = j + 1;
117 g = (d[l] - p) / (2.0 * e[l-1]);
118 r = sqrt(g * g + 1.0);
119 g = d[m-1] - p + e[l-1] / (g + fabs(r) * sign(g));
120 s = 1.0;
121 c = 1.0;
122 p = 0.0;
123 mml = m - l;
124
125 for (ii = 1; ii <= mml; ii++)
126 {
127 i = m - ii;
128 f = s * e[i-1];
129 b = c * e[i-1];
130
131 if (fabs(g) <= fabs(f)) {
132 c = g / f;
133 r = sqrt(c * c + 1.0);
134 e[i] = f * r;
135 s = 1.0 / r;
136 c = c * s;
137 } else {
138 s = f / g;
139 r = sqrt(s * s + 1.0);
140 e[i] = g * r;
141 c = 1.0 / r;
142 s = s * c;
143 }
144 g = d[i] - p;
145 r = (d[i-1] - g) * s + 2.0 * c * b;
146 p = s * r;
147 d[i] = g + p;
148 g = c * r - b;
149 f = z[i];
150 z[i] = s * z[i-1] + c * f;
151 z[i-1] = c * z[i-1] - s * f;
152 }
153 d[l-1] = d[l-1] - p;
154 e[l-1] = g;
155 e[m-1] = 0.0;
156 }
157 }
158
159 // Sorting.
160 for (ii = 2; ii <= m; ii++) {
161 i = ii - 1;
162 k = i;
163 p = d[i-1];
164
165 for (j = ii; j <= n; j++) {
166 if (d[j-1] < p) {
167 k = j;
168 p = d[j-1];
169 }
170 }
171
172 if (k != i) {
173 d[k-1] = d[i-1];
174 d[i-1] = p;
175 p = z[i-1];
176 z[i-1] = z[k-1];
177 z[k-1] = p;
178 }
179 }
180}
181
182
183 void gaussLaguerre(size_t n, std::vector<double>& abscissae, DataVector<double>& weights, double scale)
184{
185 std::unique_ptr<double[]> work(new double[n]);
186
187 abscissae.resize(n);
188 weights.reset(n);
189
190 for (size_t i = 0; i != n; ++i) {
191 abscissae[i] = double(2 * i + 1);
192 work[i] = double(i + 1);
193 }
194
195 std::fill(weights.begin(), weights.end(), 0.);
196 weights[0] = 1.0;
197
198 imtqlx(int(n), &abscissae.front(), work.get(), weights.data());
199
200 double iscale = 1. / scale;
201 for (size_t i = 0; i < n; ++i) {
202 double factor = exp(abscissae[i]);
203 if (!std::isfinite(factor)) factor = 0.;
204 weights[i] *= iscale * weights[i] * factor;
205 abscissae[i] *= iscale;
206 }
207}
208
209}}} // # namespace plask::optical::modal