#include <cmath>
#include <iostream>
#include "nr3.h"
#include "quadrature.h"
#include "interp_1d.h"
#include "romberg.h"
#include "gamma.h"
#include "gauss_wgts.h"

/*
 * define some parameters that are used in separate routines.
 * Since the functions called in the integration algorithms can
 * take only  1 argument, we have to specify the other arguments
 * globally.
*/

long int funccalls; // to check how many calculations are needed.

/*
 * try different integration modes:
 * -1 Gaussian quadrature, integral over z done analytically
 * 0,1 Gaussian quadrature, 40 or 50 steps in z
 * 2 Simpsons rule
 * 3 Romberg integration
*/
int mode; 

/*
 * integrand = 2  - integrate over dr from 0 to 1.6
 * integrand = 1 -  intergrate over dtheta from 0 to pi/2 (dr/dtheta = 4.8 sin^theta cos theta
 * integrand =0 - integrate over dx = e^theta from -80 to ln(PI/2), dx/dtheta = theta
 *
*/
int integrand; 

const double PI = acos(-1); // pi to machine accuracy

double global_r; // we pass the value r to the function that calculated the inner integral over z 

int PREC1a = 40; // number of Gaussian-quadrature abscissa for the integration over z

int PREC1b = 50; // number of Gaussian-quadrature abscissa for the integration over z

double EPS=1e-9; // accuracy outer integral

VecDoub x1a(PREC1a),weight1a(PREC1a); // for Gaussian quadrature 
VecDoub x1b(PREC1b),weight1b(PREC1b); // for Gaussian quadrature

//---------------------------------------------------------------------//

/*
 * functions for the integration routines
 * func3 returns r times f(r,z). z is passed as argument
*/
Doub func3(Doub z) {
   funccalls++;
   if (global_r<1e-100) global_r = 1e-100; // z can be smaller than -1.69
   return pow(global_r,z+1);
}

/*
 * func_r returns the function used to take a step in the integral over dr
*/

Doub func_r(Doub r) {
   if (funccalls>2e10) throw; // less then 20 billion trials, to keep this manageable.
   global_r = r;
   double t=pow((r/1.6),1/3.0);
   if (t>1.0) t=1.0; // catch round-off errors
   t=asin(t);
   double result = 2 * PI;
   double zmax = 1.3*cos(t)-0.5*cos(2*t)-0.2*cos(3*t)-0.1*cos(4*t);
   double zmin = -1.3*cos(t)-0.5*cos(2*t)+0.2*cos(3*t)-0.1*cos(4*t);
   if (zmax<zmin) cerr << " check! zmax < zmin for t = " << t << endl;
   if (zmax-zmin<1e-10) return (zmax-zmin)*func3(zmax);

   if (mode==-1) {
	funccalls++;
 	return result*(pow(r,zmax+1)-pow(r,zmin+1))/log(r);
   }

   if (mode==0) { // Gaussian quadrature, PREC1a steps
// use z = zmin + (zmax-zmin) x1, x1 runs from 0 to 1
// in this manner, you only have to calculate the weights with Gauleg once
         double sum=0;
         for (int i=0; i<PREC1a; i++) {
	   double z=zmin+x1a[i]*(zmax-zmin);
           sum+=weight1a[i]*func3(z);
         }
         return result*sum*(zmax-zmin);
    }
    if (mode==1) { // Gaussian quadrature, PREC1b steps
         double sum=0;
         for (int i=0; i<PREC1b; i++) {
	   double z=zmin+x1b[i]*(zmax-zmin);
           sum+=weight1b[i]*func3(z);
         }
         return result*sum*(zmax-zmin);
    }
    if (mode==2) return result*qsimp(func3,zmin,zmax,EPS);
    return result*qromb(func3,zmin,zmax,EPS);
}

Doub func_t(Doub t) {
	double r=1.6*pow(sin(t),3);	
	double drdt=sin(t)*sin(t)*cos(t)*4.8;
	return drdt*func_r(r);
}

Doub func_x(Doub x) {
	double t=exp(x);
	return t*func_t(t);
}

Doub integral(double x) {
	if (integrand==0) return func_x(x);
	if (integrand==1) return func_t(x);
	if (integrand==2) return func_r(x);
	cerr << " unforeseen integrand ";
	return 0;
}

//-----------------------------------------------------------------------//

int main() {
	gauleg(0,1,x1a,weight1a); // calculate the Gauss-legendre weights
	gauleg(0,1,x1b,weight1b); // calculate the Gauss-legendre weights
	cout << setprecision(16);
	for (integrand=2; integrand<3; integrand++) {
// try integration over dx = d e^theta, over dtheta, and over dr
		double xmin,xmax;
		if (integrand==0) {
			xmin=-80;
			xmax=log(PI/2);
			cout << "integration over dx = exp(theta), from theta=e-80 to theta=pi/2 " << endl;
		}
		if (integrand==1) {
			xmin=0;
			xmax=PI/2;
			cout << "integration over d(theta), from theta=1e-11 to theta=pi/2 " << endl;
		}
		if (integrand==2) {
			xmin=0; 
			xmax=1.6;
			cout << "integration over dr, from 1e-32 to 1.6 " << endl;
		}
        	for (mode=-1; mode<4; mode++) {
	    		if (mode==-1) cout << " integral over z done analytically " << endl;
	    		if (mode==0) cout << " Gauss-Legendre, 40 steps in z " << endl; 
	    		if (mode==1) cout << " Gauss-Legendre, 50 steps in z " << endl; 
            		if (mode<2) {
				for (int PREC2=50; PREC2<1000;PREC2+=50) {
	    				funccalls=0;
					VecDoub x2(PREC2),weight2(PREC2);
					gauleg(xmin,xmax,x2,weight2);
					double sum=0;
					for (int i=0; i<PREC2;i++) sum+=weight2[i]*integral(x2[i]);
					cout << "  " <<  PREC2<< " steps in Gaussian quadrature, integral: " <<  sum << " , " << funccalls << " function calls" << endl;
        			}
			}
			if (mode==2) {
				try {
				if (integrand==1) xmin=1e-30;
				if (integrand==2) xmin=1e-10;
				funccalls=0;
				EPS=1e-4;
				cout <<  "Simpson, EPS = " << EPS << " result ";
				cout << qsimp(integral,xmin,xmax,EPS);
				cout << " calls " << funccalls << endl;
/*
				EPS=1e-5;
				cout <<  "Simpson, EPS = " << EPS << " result ";
				cout << qsimp(integral,xmin,xmax,EPS);
				cout << " calls " << funccalls << endl;
*/
				}
				catch (...) {
					cout << " too many calls : " << funccalls << endl;
				}
			}
			if (mode==3) {
				try {
				funccalls=0;
				EPS=1e-6;
				cout <<  "Romberg, EPS = " << EPS << " result ";
				cout << qromb(integral,xmin,xmax,EPS);
				cout << " calls " << funccalls << endl;
				EPS=1e-8;
				cout <<  "Romberg, EPS = " << EPS << " result ";
				cout << qromb(integral,xmin,xmax,EPS);
				cout << " calls " << funccalls << endl;
				}
				catch (...) {
					cout << " too many calls : " << funccalls << endl;
				}
			}
		}
		cout << endl;
	}
	return 0;
}