defint.h

/********************************************************************/
/* Copyright (c) 2015 DG Technologies Inc. and Yuzi Mizuno          */
/* All rights reserved.                                             */
/********************************************************************/
#ifndef _MGDefint_HH_
#define _MGDefint_HH_

#include <math.h>
#include "mg/defintArea.h"

template<class func>
double mgDefint(
    func& f,
    double a,   
    double b,   
    double eps, 
    int l=0 
){
    if(l!=0 && l!=1)
        return 0.;

    double d;
    int i;
    double wm, wp, fac;

    MGDefintArea::instance().init();//Check if defintArea is initialized.
    int npow=MGDefintArea::instance().m_npow, nend=MGDefintArea::instance().m_nend;
    double *a0=MGDefintArea::instance().m_a0,
        *am=MGDefintArea::instance().m_am,
        *ap=MGDefintArea::instance().m_ap;
    double b0=MGDefintArea::instance().m_b0, *bb=MGDefintArea::instance().m_bb;
    double eps0=MGDefintArea::instance().m_eps0;

    double shfm, shfp;
    fac = (b-a)*.5;
    if(l == 0){
        shfm = (b+a)*.5;
        shfp = shfm;
    }else{
        shfm = 0.;
        shfp = 0.;
    }

    const int lp1 = l + 1;
    const int id=lp1*mgdefintlen-(mgdefintlen+1);
    double epsv;
    eps=fabs(eps);
    if(eps>=eps0)
        epsv = eps;
    else epsv = eps0;
    double epsq = sqrt(epsv) * .2;

    d = a0[l] * fac + shfp;
    double vnew = f(d) * b0;

/*     ---- initial step ---- */
/*          integrate with mesh size = 0.5 and check decay of integrand */
    double h = .5;
    int is=1; is=is<<npow;
    int im = is;
    int km=0, kp=0, nm=0, np=0;
   for(i=is; im<0 ? i>=nend : i<=nend; i+=im){

    if(km<=1){
        d = am[i+id]*fac + shfm;
        wm = f(d) * bb[i-1];
        vnew += wm;
        if(fabs(wm)<=epsv){
            ++km;
            if (km >= 2) nm = i - im;
        }else
            km = 0;
    }
    if (kp<=1) {
        d = ap[i+id]*fac + shfp;
        wp = f(d) * bb[i-1];
        vnew += wp;
        if(fabs(wp)<=epsv){
            ++kp;
            if (kp>=2) np=i-im;
        }else
            kp = 0;
    }
    if(km==2 && kp==2)
        break;

    }

    double epsm = 0., epsp = 0.;
    if (nm == 0) {
        nm = nend;
        epsm = sqrt(fabs(wm));
    }

    if(np == 0){
        np = nend;
        epsp = sqrt(fabs(wp));
    }

    if (epsq<epsm)
        epsq = epsm;
    if (epsq<epsp)
        epsq = epsp;

/*     ---- general step ---- */
    double vold = h * fac * vnew;
    for(int mstep = 1; mstep<=npow; ++mstep){
        vnew = 0.;
        int ih = is;
        is /= 2;

        for(i=is; ih<0 ? i>=nm : i<=nm; i+=ih) {
            d = am[i+id] * fac + shfm;
            vnew += f(d) * bb[i - 1];
        }

        for(i=is; ih<0 ? i>=np : i<=np; i+=ih) {
            d = ap[i+id] * fac + shfp;
            vnew += f(d) * bb[i-1];
        }

        vnew = (vold + h*fac*vnew) * .5;
        if((d=vnew-vold,fabs(d))<epsq)
            return vnew;//converged and return

        h *= .5;
        vold = vnew;
    }
    return vnew;
}
 // end of ALGORITHM group

#endif