PLaSK library

Loading...
Searching...
No Matches
patterson.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 "patterson.hpp"
15#include "patterson-data.hpp"
16
17 namespace plask { namespace optical { namespace modal {
18
19 template <typename S, typename T>
20S patterson(const std::function<S(T)>& fun, T a, T b, double& err, unsigned* order)
21{
22 double eps = err;
23 err *= 2.;
24
25 S result, result2;
26 T D = (b - a) / 2., Z = (a + b) / 2.;
27
28 S values[256]; std::fill_n(values, 256, 0.);
29 values[0] = fun(Z);
30 result = (b - a) * values[0];
31
32 unsigned n;
33
34 for (n = 1; err > eps && n < 9; ++n) {
35 const unsigned N = 1 << n; // number of points in current iteration on one side of zero
36 const unsigned stp = 256 >> n; // step in x array
37 result2 = result;
38 // Compute integral on scaled range [-1, +1]
39 result = patterson_weights[n][0] * values[255];
40 for (unsigned i = stp, j = 1; j < N; i += stp, ++j) {
41 if (j % 2) { // only odd points are new ones
42 const double x = patterson_points[i];
43 T z1 = Z - D * x;
44 T z2 = Z + D * x;
45 values[i] = fun(z1) + fun(z2);
46 }
47 result += patterson_weights[n][j] * values[i];
48 }
49 // Rescale integral and compute error
50 result *= D;
51 err = abs(1. - result2 / result);
52 }
53
54#ifndef NDEBUG
55 writelog(LOG_DEBUG, "Patterson quadrature for {0} points, error = {1}", (1<<n)-1, err);
56#endif
57
58 if (order) *order = n - 1;
59
60 return result;
61}
62
63 template double patterson<double,double>(const std::function<double(double)>& fun, double a, double b, double& err, unsigned* order);
64
65}}} // namespace plask::optical::effective