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 effective {
18
19 template <typename S, typename T>
20S patterson(const std::function<S(T)>& fun, T a, T b, double& err)
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];
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 unsigned N = 1 << n; // number of points in current iteration on one side of zero
36 unsigned N0 = N >> 1; // number of points in previous iteration on one side of zero
37 result2 = result;
38 result = 0.;
39 // Add previously computed points
40 for (unsigned i = 0; i < N0; ++i) {
41 result += patterson_weights[n][i] * values[i];
42 }
43 // Compute and add new points
44 for (unsigned i = N0; i < N; ++i) {
45 double x = patterson_points[i];
46 T z1 = Z - D * x;
47 T z2 = Z + D * x;
48 values[i] = fun(z1) + fun(z2);
49 result += patterson_weights[n][i] * values[i];
50 }
51 // Rescale integral and compute error
52 result *= D;
53 err = abs(1. - result2 / result);
54 }
55
56#ifndef NDEBUG
57 writelog(LOG_DEBUG, "Patterson quadrature for {0} points, error = {1}", (2<<n)-1, err);
58#endif
59
60 return result;
61}
62
63 template double patterson<double,double>(const std::function<double(double)>& fun, double a, double b, double& err);
64
65}}} // namespace plask::optical::effective