//  RombergDemo.java
//  Demonstrates Romberg integration.
//  Avoid premature exit: use MINITER>3, TOL<1E-8.
//  funcparms[] passes straight through to func(). 
//  M.Lampton UCB SSL 2009

public class RombergDemo
{
    static int ncalls = 0; 

    public static void main(String[] args)
    {
        double funcparms[] = {10.0};
        double ans = doRomberg(0.0, 1.0, funcparms);
        double err = Math.abs(1 - ans); 
        System.out.println(" Romberg err = "+err+";   ncalls = "+ncalls); 
    }

    static double func(double x, double parms[])
    // This is a negative exponential, normalized to unit integral.
    // Its parameter controls the concentration near the origin.
    {
        ncalls++; 
        double coef = parms[0]/(1 - Math.exp(-parms[0])); 
        return coef*Math.exp(-x*parms[0]);
    } 

    //-----Romberg package------

    static int MAXITER = 25; 
    static int MINITER = 3;  
    static double r[][] = new double[MAXITER][MAXITER]; 
    static double TOL = 1E-9;
    static double ERRCODE = 999.9; 

    static double doRomberg(double a, double b, double fp[])
    // http://en.wikipedia.org/wiki/Romberg%27s_method
    {
        r[0][0] = 0.5*(b-a)*(func(a, fp)+func(b, fp)); 
        int twon = 1; 
        for (int n=1; n<MAXITER; n++)
        {
            twon *= 2; 
            int kmax = twon/2; 
            double hn = (b-a)/twon; 
            r[n][0] = 0.5*r[n-1][0]; 
            for (int k=1; k<=kmax; k++)
              r[n][0] += hn * func(a+(2*k-1)*hn, fp); 
            for (int m=1; m<=n; m++)
            {
                double fourm = Math.pow(4,m); 
                r[n][m] = (fourm*r[n][m-1] - r[n-1][m-1])/(fourm-1); 
            }
            double err = Math.abs(r[n][n] - r[n][n-1]); 
            if ((n>=MINITER) && (err < TOL))
              return r[n][n]; 
        }
        return ERRCODE;
    }        
}