PLaSK library

Loading...
Searching...
No Matches
gauss_legendre.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_legendre.hpp"
15#include "fortran.hpp"
16
17#include <boost/math/special_functions/legendre.hpp>
18 using boost::math::legendre_p;
19
20 namespace plask { namespace optical { namespace modal {
21
22
23 void gaussLegendre(size_t n, std::vector<double>& abscissae, DataVector<double>& weights)
24{
25 int info;
26
27 abscissae.assign(n, 0.);
28 weights.reset(n);
29
30 for (size_t i = 1; i != n; ++i)
31 weights[i-1] = 0.5 / std::sqrt(1. - 0.25/double(i*i));
32
33 dsterf(int(n), &abscissae.front(), weights.data(), info);
34 if (info < 0) throw CriticalException("Gauss-Legendre abscissae: argument {:d} of DSTERF has bad value", -info);
35 if (info > 0) throw ComputationError("Gauss-Legendre abscissae", "could not converge in {:d}-th element", info);
36
37 double nn = double(n*n);
38 auto x = abscissae.begin();
39 auto w = weights.begin();
40 for (; x != abscissae.end(); ++x, ++w) {
41 double P = legendre_p(int(n-1), *x);
42 *w = 2. * (1. - (*x)*(*x)) / (nn * P*P);
43 }
44}
45
46}}} // # namespace plask::optical::modal