PLaSK library

Loading...
Searching...
No Matches
muller.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 "muller.hpp"
15#include "rootdigger-impl.hpp"
16#include "solver.hpp"
17 using namespace std;
18
19 namespace plask { namespace optical { namespace modal {
20
21//**************************************************************************
23dcomplex RootMuller::find(dcomplex start)
24{
25 dcomplex first = start - 0.5 * params.initial_dist;
26 dcomplex second = start + 0.5 * params.initial_dist;
27
28 writelog(LOG_DETAIL, "Searching for the root with Muller method between {0} and {1}", str(first), str(second));
29 log_value.resetCounter();
30
31 double xtol2 = params.tolx * params.tolx;
32 double fmin2 = params.tolf_min * params.tolf_min;
33 double fmax2 = params.tolf_max * params.tolf_max;
34
35 dcomplex x2 = first, x1 = second, x0 = start;
36
37 dcomplex f2 = valFunction(x2); log_value(x2, f2);
38 dcomplex f1 = valFunction(x1); log_value(x1, f1);
39 dcomplex f0 = valFunction(x0); log_value.count(x0, f0);
40
41 for (int i = 0; i < params.maxiter; ++i) {
42 if (isnan(real(f0)) || isnan(imag(f0)))
43 throw ComputationError(solver.getId(), "computed value is NaN");
44 dcomplex q = (x0 - x1) / (x1 - x2);
45 dcomplex A = q * f0 - q*(q+1.) * f1 + q*q * f2;
46 dcomplex B = (2.*q+1.) * f0 - (q+1.)*(q+1.) * f1 + q*q * f2;
47 dcomplex C = (q+1.) * f0;
48 dcomplex S = sqrt(B*B - 4.*A*C);
49 x2 = x1; f2 = f1;
50 x1 = x0; f1 = f0;
51 x0 = x1 - (x1-x2) * ( 2.*C / std::max(B+S, B-S, [](const dcomplex& a, const dcomplex& b){return abs2(a) < abs2(b);}) );
52 f0 = valFunction(x0); log_value.count(x0, f0);
53 if (abs2(f0) < fmin2 || (abs2(x0-x1) < xtol2 && abs2(f0) < fmax2)) {
54 writelog(LOG_RESULT, "Found root at " + str(x0));
55 return x0;
56 }
57 }
58
59 throw ComputationError(solver.getId(), "Muller: {0}: maximum number of iterations reached", log_value.chartName());
60}
61
62}}} // namespace plask::optical::modal