KLFunctions.hh

#ifndef KLFunctions_h
#define KLFunctions_h
 
namespace Dune{
    namespace Composites{
 
        template<class X>
            void inline evaluate_eigenValues(double ell,  X& lambda, X& freq){
                double c = 1. / ell;
                for (int i = 0; i < lambda.size(); i++){
                    lambda[i] = 2.0 * c /(freq[i] * freq[i] + c * c);
                }
            }
 
 
        double res(double x, double ell){
            double c = 1.0 / ell;
            double g = std::tan(x) - 2.0 * c * x / (x * x - c * c);
            return g;
        }
 
 
        double findRoot(double ell, double a, double b){
            double fa, fb, x, fx, error, m;
            error = 1.0;
            int iter = 0;
            while (error > 1e-6){
                fa = res(a,ell);
                fb = res(b,ell);
                m = (fb - fa) / (b - a);
                x = a - fa / m;
                fx = res(x,ell);
                if (((fa < 0) & (fx < 0)) | ((fa > 0) & (fx > 0))) { a = x; }
                else { b = x; }
                error = std::abs(fx);
                if(iter > 100000){
                    std::cout << "Warning: root finder hasn't converged. Check number of KL modes and correlation length." << std::endl;
                    break;
                }
                iter++;
            }
            return x;
        }
 
        void rootFinder(int M, double ell, std::vector<double>& answer){
            double c = 1.0 / ell;
            std::vector<double> freq(M+2);
 
            // For all intervals
 
            int m = -1;
            //Should go to M because first interval is added additionally!
            for (int i = 0; i < M+1; i++){
 
                //  std::cout << "Root i = " << i << std::endl;
 
                double w_min = (i - 0.4999) * M_PI;
                double w_max = (i + 0.4999) * M_PI;
 
                //  std::cout << "w_min = " << w_min << std::endl;
                //std::cout << "w_max = " << w_max << std::endl;
                if ((w_min <= c) && (w_max >= c)){
 
                    // If not first interval look for solution near left boundary
                    if (w_min > 0.0){
                        m += 1;
                        freq[m] = findRoot(ell,w_min,0.5*(c+w_min));
                    }
                    // Always look for solution near right boundary
                    m += 1;
                    freq[m] = findRoot(ell,0.5*(c + w_max),w_max);
                }
                else{
                    m += 1;
                    freq[m] = findRoot(ell,w_min,w_max);
                }
 
            }
 
            for (int i = 0; i < M; i++){
                answer[i] = freq[i+1];
            }
 
        }
 
    }
}
 
#endif /* KLFunctions_h */