/*
 Integration of a function using following approximations:
 # Trapezoid approximation
 # Simpson's rule
 # Monte Carlo method
 # 8-panel newton-cotes rule
 time of integration is also calculated and printed
*/
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include "trapezoid.cpp"
#include "simpson.cpp"
#include "quanc8.cpp"
#include "mc_int1.cpp"
using namespace std;

//double mc_int1(double(*)(double), double, double, int);
double f(double); 

int main()
{
    double a, b;
    double trap, simp, mc;
    double nc8, errest, flag;
    int i, n, nofun;
    int ntimes;
    const double pi = 3.1415926;
    const double abserr=0.0, relerr=1.0e-12;  //see quanc8.cpp
    double tail, errmc;    
       
    cout.precision(6);
    cout.setf(ios::fixed | ios::showpoint);

// current local time
    time_t rawtime;
    time ( &rawtime );
    cout << " The current local time is:" << endl <<" "<< ctime (&rawtime) << endl;

// current time in seconds (begin calculations)
    time_t seconds_i;
    seconds_i = time (NULL);

/* initial information */    
    a = 0.0;              // left endpoint
    b = pi;               // right endpoint
    n = 2;                // initial number of intervals
    ntimes = 16;          // number of interval doublings with nmax=2^ntimes

/* "tail part for the 0->infinity integral */
//    tail = 4.0/sqrt(b);  // it works ONLY for the last integral in the list
    tail = 0.0;


    cout << " Intervals  "<<"Trapez.   "<<"Simpson   "<< "MC      " 
         << "  MC error " << endl;

/* step 2: integration with trapezoid and simpson */
    for (i=0; i <=ntimes; i=i+1)
    {
        trap = trapezoid(f, a, b, n);
        simp =   simpson(f, a, b, n);
        mc =     mc_int1(f, errmc, a, b, n);
        cout << setw(10)  << n  << setw(10) << trap +tail  << setw(10) << simp+tail 
             << setw(10)  << mc << setw(10) << errmc << endl;
        n = n*2;
    }

/* step 3: 8-panel newton-cotes rule */
    quanc8(f, a, b, abserr, relerr, nc8, errest, nofun, flag);

    cout << endl     << "    result from quanc8 " << endl;
    cout << setw(12) << "nofun = "    << setw(10) << nofun << endl;
    cout << setw(12) << "integral = " << setw(10) << nc8+tail   << endl;
//    cout << setw(12) << "est.error = "<< setw(10) << errest << endl;
//    cout << setw(12) << "flag = "     << setw(10) << flag   << endl;    

// current time in seconds (end of calculations)
    time_t seconds_f;
    seconds_f = time (NULL);
    cout << endl << "total elapsed time = " << seconds_f - seconds_i << " seconds" << endl << endl;
  
    return 0;
}
/*
 *  Function for integration
*/
    double f(double x)
{
    const double pi = 3.1415926;
    double y;
//    y = sin(x);
//    y = cos(x);
    y = 2.0*sqrt(x)/((x+1.0)*(x+1.0));
//    y = 0.2/(pow((x-6.0),2.0) + 0.02);
//    y = x*cos(10.0*pow(x,2.0))/(pow(x,2.0) + 1.0);
//    y = sqrt(x)/(x*x + 1.0);
//    y = sin(pi*x*x)/((x-pi)*(x-pi)+1);
    return y;
}