#ifndef NUMERICAL_INTEGRATION_H
#define NUMERICAL_INTEGRATION_H
double Integrate2DLaguerre32Legendre32(double (*func)(double,double), double ay, double by)
{
  // Integrate 2D function from 0 to infinity using Gauss-Laguerre integration
  // with 32 points and from ay to by using Gauss-Legendre integration with 32 points
  //

  double xleg[32];
  double wleg[32];
  double x[32];
  double w[32];

  xleg[0] = -0.997263861849; 
  xleg[1] = -0.985611511545;
  xleg[2] = -0.964762255588; 
  xleg[3] = -0.934906075938; 
  xleg[4] = -0.896321155766; 
  xleg[5] = -0.849367613733; 
  xleg[6] = -0.794483795968; 
  xleg[7] = -0.732182118740; 
  xleg[8] = -0.663044266930; 
  xleg[9] = -0.587715757241; 
  xleg[10] = -0.506899908932; 
  xleg[11] = -0.421351276131; 
  xleg[12] = -0.331868602282; 
  xleg[13] = -0.239287362252; 
  xleg[14] = -0.144471961583; 
  xleg[15] = -0.048307665688; 
  xleg[16] = 0.048307665688; 	
  xleg[17] = 0.144471961583; 	
  xleg[18] = 0.239287362252; 	
  xleg[19] = 0.331868602282; 	
  xleg[20] = 0.421351276131; 	
  xleg[21] = 0.506899908932; 	
  xleg[22] = 0.587715757241;	
  xleg[23] = 0.663044266930; 	
  xleg[24] = 0.732182118740; 	
  xleg[25] = 0.794483795968; 	
  xleg[26] = 0.849367613733; 	
  xleg[27] = 0.896321155766; 	
  xleg[28] = 0.934906075938; 	
  xleg[29] = 0.964762255588; 	
  xleg[30] = 0.985611511545; 	
  xleg[31] = 0.997263861849; 	

  wleg[0] = 0.007018610009;
  wleg[1] = 0.016274394716;
  wleg[2] = 0.025392065309;
  wleg[3] = 0.034273862913;
  wleg[4] = 0.042835898022;
  wleg[5] = 0.050998059262;
  wleg[6] = 0.058684093479;
  wleg[7] = 0.065822222776;
  wleg[8] = 0.072345794109;
  wleg[9] = 0.078193895787;
  wleg[10] = 0.083311924227;
  wleg[11] = 0.087652093004;
  wleg[12] = 0.091173878696;
  wleg[13] = 0.093844399081;
  wleg[14] = 0.095638720079;
  wleg[15] = 0.096540088515;
  wleg[16] = 0.096540088515;
  wleg[17] = 0.095638720079;
  wleg[18] = 0.093844399081;
  wleg[19] = 0.091173878696;
  wleg[20] = 0.087652093004;
  wleg[21] = 0.083311924227;
  wleg[22] = 0.078193895787;
  wleg[23] = 0.072345794109;
  wleg[24] = 0.065822222776;
  wleg[25] = 0.058684093479;
  wleg[26] = 0.050998059262;
  wleg[27] = 0.042835898022;
  wleg[28] = 0.034273862913;
  wleg[29] = 0.025392065309;
  wleg[30] = 0.016274394716;
  wleg[31] = 0.007018610009;

  double xlag[32];
  double wlag[32];

  xlag[0] = 0.044489365833; 	  	
  xlag[1] = 0.234526109520;	 	
  xlag[2] = 0.576884629302; 		
  xlag[3] = 1.072448753818; 		
  xlag[4] = 1.722408776445; 		
  xlag[5] = 2.528336706426; 		
  xlag[6] = 3.492213273022; 		
  xlag[7] = 4.616456769750; 		
  xlag[8] = 5.903958504174; 		
  xlag[9] = 7.358126733186; 		
  xlag[10] = 8.982940924213; 		
  xlag[11] = 10.783018632540; 	
  xlag[12] = 12.763697986743; 	
  xlag[13] = 14.931139755523; 	
  xlag[14] = 17.292454336715; 	
  xlag[15] = 19.855860940336; 	
  xlag[16] = 22.630889013197; 	
  xlag[17] = 25.628636022459; 	
  xlag[18] = 28.862101816323; 	
  xlag[19] = 32.346629153965; 	
  xlag[20] = 36.100494805752; 	
  xlag[21] = 40.145719771539; 	
  xlag[22] = 44.509207995755; 	
  xlag[23] = 49.224394987309; 	
  xlag[24] = 54.333721333397; 	
  xlag[25] = 59.892509162134; 	
  xlag[26] = 65.975377287935; 	
  xlag[27] = 72.687628090663;	
  xlag[28] = 80.187446977914;	
  xlag[29] = 88.735340417892; 	
  xlag[30] = 98.829542868284;	
  xlag[31] = 111.751398097938; 	

  wlag[0] = 0.114187105768;
  wlag[1] = 0.266065216898;
  wlag[2] = 0.418793137325;
  wlag[3] = 0.572532846500;
  wlag[4] = 0.727648788381;
  wlag[5] = 0.884536719340;
  wlag[6] = 1.043618875892;
  wlag[7] = 1.205349274152;
  wlag[8] = 1.370221338522;
  wlag[9] = 1.538777256469;
  wlag[10] = 1.711619352686;
  wlag[11] = 1.889424063449;
  wlag[12] = 2.072959340247;
  wlag[13] = 2.263106633997;
  wlag[14] = 2.460889072488;
  wlag[15] = 2.667508126397;
  wlag[16] = 2.884392092922;
  wlag[17] = 3.113261327040;
  wlag[18] = 3.356217692596;
  wlag[19] = 3.615869856484;
  wlag[20] = 3.895513044949;
  wlag[21] = 4.199394104712;
  wlag[22] = 4.533114978534;
  wlag[23] = 4.904270287611;
  wlag[24] = 5.323500972024;
  wlag[25] = 5.806333214234;
  wlag[26] = 6.376614674160;
  wlag[27] = 7.073526580707;
  wlag[28] = 7.967693509296;
  wlag[29] = 9.205040331278;
  wlag[30] = 11.163013090768;
  wlag[31] = 15.390180415261;

  double sum = 0.;

  for(int i = 0 ; i < 32 ; i++){
    for(int j = 0 ; j < 32 ; j++){
      x[j] = (by-ay)/2.*xleg[j] + (by+ay)/2.;
      w[j] = (by-ay)/2.*wleg[j];

      sum += wlag[i]*w[j]*func(xlag[i],x[j]);
    }
  }
  return sum;

}

void GetCoefs2DLaguerre32Legendre32(double ay, double by, std::vector<double> &xlag, std::vector<double> &wlag, std::vector<double> &xleg, std::vector<double> &wleg) {
	  xlag.resize(32);
	  wlag.resize(32);
	  xleg.resize(32);
	  wleg.resize(32);
	
	  xleg[0] = -0.997263861849; 
	  xleg[1] = -0.985611511545;
	  xleg[2] = -0.964762255588; 
	  xleg[3] = -0.934906075938; 
	  xleg[4] = -0.896321155766; 
	  xleg[5] = -0.849367613733; 
	  xleg[6] = -0.794483795968; 
	  xleg[7] = -0.732182118740; 
	  xleg[8] = -0.663044266930; 
	  xleg[9] = -0.587715757241; 
	  xleg[10] = -0.506899908932; 
	  xleg[11] = -0.421351276131; 
	  xleg[12] = -0.331868602282; 
	  xleg[13] = -0.239287362252; 
	  xleg[14] = -0.144471961583; 
	  xleg[15] = -0.048307665688; 
	  xleg[16] = 0.048307665688; 	
	  xleg[17] = 0.144471961583; 	
	  xleg[18] = 0.239287362252; 	
	  xleg[19] = 0.331868602282; 	
	  xleg[20] = 0.421351276131; 	
	  xleg[21] = 0.506899908932; 	
	  xleg[22] = 0.587715757241;	
	  xleg[23] = 0.663044266930; 	
	  xleg[24] = 0.732182118740; 	
	  xleg[25] = 0.794483795968; 	
	  xleg[26] = 0.849367613733; 	
	  xleg[27] = 0.896321155766; 	
	  xleg[28] = 0.934906075938; 	
	  xleg[29] = 0.964762255588; 	
	  xleg[30] = 0.985611511545; 	
	  xleg[31] = 0.997263861849; 	

	  wleg[0] = 0.007018610009;
	  wleg[1] = 0.016274394716;
	  wleg[2] = 0.025392065309;
	  wleg[3] = 0.034273862913;
	  wleg[4] = 0.042835898022;
	  wleg[5] = 0.050998059262;
	  wleg[6] = 0.058684093479;
	  wleg[7] = 0.065822222776;
	  wleg[8] = 0.072345794109;
	  wleg[9] = 0.078193895787;
	  wleg[10] = 0.083311924227;
	  wleg[11] = 0.087652093004;
	  wleg[12] = 0.091173878696;
	  wleg[13] = 0.093844399081;
	  wleg[14] = 0.095638720079;
	  wleg[15] = 0.096540088515;
	  wleg[16] = 0.096540088515;
	  wleg[17] = 0.095638720079;
	  wleg[18] = 0.093844399081;
	  wleg[19] = 0.091173878696;
	  wleg[20] = 0.087652093004;
	  wleg[21] = 0.083311924227;
	  wleg[22] = 0.078193895787;
	  wleg[23] = 0.072345794109;
	  wleg[24] = 0.065822222776;
	  wleg[25] = 0.058684093479;
	  wleg[26] = 0.050998059262;
	  wleg[27] = 0.042835898022;
	  wleg[28] = 0.034273862913;
	  wleg[29] = 0.025392065309;
	  wleg[30] = 0.016274394716;
	  wleg[31] = 0.007018610009;

	  xlag[0] = 0.044489365833; 	  	
	  xlag[1] = 0.234526109520;	 	
	  xlag[2] = 0.576884629302; 		
	  xlag[3] = 1.072448753818; 		
	  xlag[4] = 1.722408776445; 		
	  xlag[5] = 2.528336706426; 		
	  xlag[6] = 3.492213273022; 		
	  xlag[7] = 4.616456769750; 		
	  xlag[8] = 5.903958504174; 		
	  xlag[9] = 7.358126733186; 		
	  xlag[10] = 8.982940924213; 		
	  xlag[11] = 10.783018632540; 	
	  xlag[12] = 12.763697986743; 	
	  xlag[13] = 14.931139755523; 	
	  xlag[14] = 17.292454336715; 	
	  xlag[15] = 19.855860940336; 	
	  xlag[16] = 22.630889013197; 	
	  xlag[17] = 25.628636022459; 	
	  xlag[18] = 28.862101816323; 	
	  xlag[19] = 32.346629153965; 	
	  xlag[20] = 36.100494805752; 	
	  xlag[21] = 40.145719771539; 	
	  xlag[22] = 44.509207995755; 	
	  xlag[23] = 49.224394987309; 	
	  xlag[24] = 54.333721333397; 	
	  xlag[25] = 59.892509162134; 	
	  xlag[26] = 65.975377287935; 	
	  xlag[27] = 72.687628090663;	
	  xlag[28] = 80.187446977914;	
	  xlag[29] = 88.735340417892; 	
	  xlag[30] = 98.829542868284;	
	  xlag[31] = 111.751398097938; 	

	  wlag[0] = 0.114187105768;
	  wlag[1] = 0.266065216898;
	  wlag[2] = 0.418793137325;
	  wlag[3] = 0.572532846500;
	  wlag[4] = 0.727648788381;
	  wlag[5] = 0.884536719340;
	  wlag[6] = 1.043618875892;
	  wlag[7] = 1.205349274152;
	  wlag[8] = 1.370221338522;
	  wlag[9] = 1.538777256469;
	  wlag[10] = 1.711619352686;
	  wlag[11] = 1.889424063449;
	  wlag[12] = 2.072959340247;
	  wlag[13] = 2.263106633997;
	  wlag[14] = 2.460889072488;
	  wlag[15] = 2.667508126397;
	  wlag[16] = 2.884392092922;
	  wlag[17] = 3.113261327040;
	  wlag[18] = 3.356217692596;
	  wlag[19] = 3.615869856484;
	  wlag[20] = 3.895513044949;
	  wlag[21] = 4.199394104712;
	  wlag[22] = 4.533114978534;
	  wlag[23] = 4.904270287611;
	  wlag[24] = 5.323500972024;
	  wlag[25] = 5.806333214234;
	  wlag[26] = 6.376614674160;
	  wlag[27] = 7.073526580707;
	  wlag[28] = 7.967693509296;
	  wlag[29] = 9.205040331278;
	  wlag[30] = 11.163013090768;
	  wlag[31] = 15.390180415261;

      //double sum = 0.;

      //for(int i = 0 ; i < 32 ; i++){
        for(int j = 0 ; j < 32 ; j++){
		  xleg[j] = (by-ay)/2.*xleg[j] + (by+ay)/2.;
		  wleg[j] = (by-ay)/2.*wleg[j];

		  //sum += wlag[i]*w[j]*func(xlag[i],x[j]);
		}
	  //}
}

void GetCoefs2DLegendre32Legendre32(double ay, double by, double a2y, double b2y, std::vector<double> &xlag, std::vector<double> &wlag, std::vector<double> &xleg, std::vector<double> &wleg) {
	  xlag.resize(32);
	  wlag.resize(32);
	  xleg.resize(32);
	  wleg.resize(32);
	
	  xleg[0] = -0.997263861849; 
	  xleg[1] = -0.985611511545;
	  xleg[2] = -0.964762255588; 
	  xleg[3] = -0.934906075938; 
	  xleg[4] = -0.896321155766; 
	  xleg[5] = -0.849367613733; 
	  xleg[6] = -0.794483795968; 
	  xleg[7] = -0.732182118740; 
	  xleg[8] = -0.663044266930; 
	  xleg[9] = -0.587715757241; 
	  xleg[10] = -0.506899908932; 
	  xleg[11] = -0.421351276131; 
	  xleg[12] = -0.331868602282; 
	  xleg[13] = -0.239287362252; 
	  xleg[14] = -0.144471961583; 
	  xleg[15] = -0.048307665688; 
	  xleg[16] = 0.048307665688; 	
	  xleg[17] = 0.144471961583; 	
	  xleg[18] = 0.239287362252; 	
	  xleg[19] = 0.331868602282; 	
	  xleg[20] = 0.421351276131; 	
	  xleg[21] = 0.506899908932; 	
	  xleg[22] = 0.587715757241;	
	  xleg[23] = 0.663044266930; 	
	  xleg[24] = 0.732182118740; 	
	  xleg[25] = 0.794483795968; 	
	  xleg[26] = 0.849367613733; 	
	  xleg[27] = 0.896321155766; 	
	  xleg[28] = 0.934906075938; 	
	  xleg[29] = 0.964762255588; 	
	  xleg[30] = 0.985611511545; 	
	  xleg[31] = 0.997263861849; 	

	  wleg[0] = 0.007018610009;
	  wleg[1] = 0.016274394716;
	  wleg[2] = 0.025392065309;
	  wleg[3] = 0.034273862913;
	  wleg[4] = 0.042835898022;
	  wleg[5] = 0.050998059262;
	  wleg[6] = 0.058684093479;
	  wleg[7] = 0.065822222776;
	  wleg[8] = 0.072345794109;
	  wleg[9] = 0.078193895787;
	  wleg[10] = 0.083311924227;
	  wleg[11] = 0.087652093004;
	  wleg[12] = 0.091173878696;
	  wleg[13] = 0.093844399081;
	  wleg[14] = 0.095638720079;
	  wleg[15] = 0.096540088515;
	  wleg[16] = 0.096540088515;
	  wleg[17] = 0.095638720079;
	  wleg[18] = 0.093844399081;
	  wleg[19] = 0.091173878696;
	  wleg[20] = 0.087652093004;
	  wleg[21] = 0.083311924227;
	  wleg[22] = 0.078193895787;
	  wleg[23] = 0.072345794109;
	  wleg[24] = 0.065822222776;
	  wleg[25] = 0.058684093479;
	  wleg[26] = 0.050998059262;
	  wleg[27] = 0.042835898022;
	  wleg[28] = 0.034273862913;
	  wleg[29] = 0.025392065309;
	  wleg[30] = 0.016274394716;
	  wleg[31] = 0.007018610009;

	  xlag[0] = 0.044489365833; 	  	
	  xlag[1] = 0.234526109520;	 	
	  xlag[2] = 0.576884629302; 		
	  xlag[3] = 1.072448753818; 		
	  xlag[4] = 1.722408776445; 		
	  xlag[5] = 2.528336706426; 		
	  xlag[6] = 3.492213273022; 		
	  xlag[7] = 4.616456769750; 		
	  xlag[8] = 5.903958504174; 		
	  xlag[9] = 7.358126733186; 		
	  xlag[10] = 8.982940924213; 		
	  xlag[11] = 10.783018632540; 	
	  xlag[12] = 12.763697986743; 	
	  xlag[13] = 14.931139755523; 	
	  xlag[14] = 17.292454336715; 	
	  xlag[15] = 19.855860940336; 	
	  xlag[16] = 22.630889013197; 	
	  xlag[17] = 25.628636022459; 	
	  xlag[18] = 28.862101816323; 	
	  xlag[19] = 32.346629153965; 	
	  xlag[20] = 36.100494805752; 	
	  xlag[21] = 40.145719771539; 	
	  xlag[22] = 44.509207995755; 	
	  xlag[23] = 49.224394987309; 	
	  xlag[24] = 54.333721333397; 	
	  xlag[25] = 59.892509162134; 	
	  xlag[26] = 65.975377287935; 	
	  xlag[27] = 72.687628090663;	
	  xlag[28] = 80.187446977914;	
	  xlag[29] = 88.735340417892; 	
	  xlag[30] = 98.829542868284;	
	  xlag[31] = 111.751398097938; 	

	  wlag[0] = 0.114187105768;
	  wlag[1] = 0.266065216898;
	  wlag[2] = 0.418793137325;
	  wlag[3] = 0.572532846500;
	  wlag[4] = 0.727648788381;
	  wlag[5] = 0.884536719340;
	  wlag[6] = 1.043618875892;
	  wlag[7] = 1.205349274152;
	  wlag[8] = 1.370221338522;
	  wlag[9] = 1.538777256469;
	  wlag[10] = 1.711619352686;
	  wlag[11] = 1.889424063449;
	  wlag[12] = 2.072959340247;
	  wlag[13] = 2.263106633997;
	  wlag[14] = 2.460889072488;
	  wlag[15] = 2.667508126397;
	  wlag[16] = 2.884392092922;
	  wlag[17] = 3.113261327040;
	  wlag[18] = 3.356217692596;
	  wlag[19] = 3.615869856484;
	  wlag[20] = 3.895513044949;
	  wlag[21] = 4.199394104712;
	  wlag[22] = 4.533114978534;
	  wlag[23] = 4.904270287611;
	  wlag[24] = 5.323500972024;
	  wlag[25] = 5.806333214234;
	  wlag[26] = 6.376614674160;
	  wlag[27] = 7.073526580707;
	  wlag[28] = 7.967693509296;
	  wlag[29] = 9.205040331278;
	  wlag[30] = 11.163013090768;
	  wlag[31] = 15.390180415261;

      //double sum = 0.;

      //for(int i = 0 ; i < 32 ; i++){
        for(int j = 0 ; j < 32 ; j++){
		  xlag[j] = (by-ay)/2.*xleg[j] + (by+ay)/2.;
		  wlag[j] = (by-ay)/2.*wleg[j];
		  xleg[j] = (b2y-a2y)/2.*xleg[j] + (b2y+a2y)/2.;
		  wleg[j] = (b2y-a2y)/2.*wleg[j];

		  //sum += wlag[i]*w[j]*func(xlag[i],x[j]);
		}
	  //}
}

void GetCoefsIntegrateLegendre32(double a, double b, std::vector<double> &x, std::vector<double> &w)
{
  // Integrate function from a to b using Legendre-Gaussian integration 
  // with 32 points.
  //

  x.resize(32);
  w.resize(32);

  x[0] = -0.997263861849; 
  x[1] = -0.985611511545;
  x[2] = -0.964762255588; 
  x[3] = -0.934906075938; 
  x[4] = -0.896321155766; 
  x[5] = -0.849367613733; 
  x[6] = -0.794483795968; 
  x[7] = -0.732182118740; 
  x[8] = -0.663044266930; 
  x[9] = -0.587715757241; 
  x[10] = -0.506899908932; 
  x[11] = -0.421351276131; 
  x[12] = -0.331868602282; 
  x[13] = -0.239287362252; 
  x[14] = -0.144471961583; 
  x[15] = -0.048307665688; 
  x[16] = 0.048307665688; 	
  x[17] = 0.144471961583; 	
  x[18] = 0.239287362252; 	
  x[19] = 0.331868602282; 	
  x[20] = 0.421351276131; 	
  x[21] = 0.506899908932; 	
  x[22] = 0.587715757241;	
  x[23] = 0.663044266930; 	
  x[24] = 0.732182118740; 	
  x[25] = 0.794483795968; 	
  x[26] = 0.849367613733; 	
  x[27] = 0.896321155766; 	
  x[28] = 0.934906075938; 	
  x[29] = 0.964762255588; 	
  x[30] = 0.985611511545; 	
  x[31] = 0.997263861849; 	

  w[0] = 0.007018610009;
  w[1] = 0.016274394716;
  w[2] = 0.025392065309;
  w[3] = 0.034273862913;
  w[4] = 0.042835898022;
  w[5] = 0.050998059262;
  w[6] = 0.058684093479;
  w[7] = 0.065822222776;
  w[8] = 0.072345794109;
  w[9] = 0.078193895787;
  w[10] = 0.083311924227;
  w[11] = 0.087652093004;
  w[12] = 0.091173878696;
  w[13] = 0.093844399081;
  w[14] = 0.095638720079;
  w[15] = 0.096540088515;
  w[16] = 0.096540088515;
  w[17] = 0.095638720079;
  w[18] = 0.093844399081;
  w[19] = 0.091173878696;
  w[20] = 0.087652093004;
  w[21] = 0.083311924227;
  w[22] = 0.078193895787;
  w[23] = 0.072345794109;
  w[24] = 0.065822222776;
  w[25] = 0.058684093479;
  w[26] = 0.050998059262;
  w[27] = 0.042835898022;
  w[28] = 0.034273862913;
  w[29] = 0.025392065309;
  w[30] = 0.016274394716;
  w[31] = 0.007018610009;
  
  for(int i = 0 ; i < 32 ; i++){
    w[i] = (b-a)/2.*w[i];
    x[i] = (b-a)/2.*x[i] + (b+a)/2.;
  }

  /*double sum = 0.;

  for(int i = 0 ; i < 32 ; i++){
    w[i] = (b-a)/2.*w[i];
    x[i] = (b-a)/2.*x[i] + (b+a)/2.;
    sum += w[i]*func->Eval(x[i]);
  }

  return sum;*/

}

void GetCoefsIntegrateLegendre10(double a, double b, std::vector<double> &x, std::vector<double> &w)
{
  // Integrate function from a to b using Legendre-Gaussian integration 
  // with 32 points.
  //

  x.resize(10);
  w.resize(10);

  x[0] = -0.973906529; 
  x[1] = -0.8650633667;
  x[2] = -0.6794095683; 
  x[3] = -0.4333953941; 
  x[4] = -0.148874339; 
  x[5] = 0.148874339; 
  x[6] = 0.433395394;
  x[7] = 0.6794095683; 
  x[8] = 0.8650633667; 
  x[9] = 0.973906529;

  w[0] = 0.0666713443;
  w[1] = 0.1494513491;
  w[2] = 0.219086363;
  w[3] = 0.2692667193;
  w[4] = 0.2955242247;
  w[5] = 0.295524225;
  w[6] = 0.2692667193;
  w[7] = 0.219086363;
  w[8] = 0.149451349;
  w[9] = 0.06667134431;
  
  for(int i = 0 ; i < 10 ; i++){
    w[i] = (b-a)/2.*w[i];
    x[i] = (b-a)/2.*x[i] + (b+a)/2.;
  }

  /*double sum = 0.;

  for(int i = 0 ; i < 32 ; i++){
    w[i] = (b-a)/2.*w[i];
    x[i] = (b-a)/2.*x[i] + (b+a)/2.;
    sum += w[i]*func->Eval(x[i]);
  }

  return sum;*/

}

void GetCoefsIntegrateLegendre5(double a, double b, std::vector<double> &x, std::vector<double> &w)
{
  // Integrate function from a to b using Legendre-Gaussian integration 
  // with 32 points.
  //

  x.resize(5);
  w.resize(5);

  x[0] = -0.906179846; 
  x[1] = -0.5384693101;
  x[2] = 0; 
  x[3] = 0.5384693101; 
  x[4] = 0.9061798459; 

  w[0] = 0.2369268851;
  w[1] = 0.4786286705;
  w[2] = 0.568888889;
  w[3] = 0.47862867;
  w[4] = 0.2369268851;
  
  for(int i = 0 ; i < 5 ; i++){
    w[i] = (b-a)/2.*w[i];
    x[i] = (b-a)/2.*x[i] + (b+a)/2.;
  }

  /*double sum = 0.;

  for(int i = 0 ; i < 32 ; i++){
    w[i] = (b-a)/2.*w[i];
    x[i] = (b-a)/2.*x[i] + (b+a)/2.;
    sum += w[i]*func->Eval(x[i]);
  }

  return sum;*/

}

void GetCoefsIntegrateLegendre40(double a, double b, std::vector<double> &x, std::vector<double> &w)
{
  // Integrate function from a to b using Legendre-Gaussian integration 
  // with 40 points.
  //

  x.resize(40);
  w.resize(40);

  x[0] = -0.998237709711; 	 
  x[1] = -0.990726238699;	 
  x[2] = -0.977259949984; 	 
  x[3] = -0.957916819214; 	 
  x[4] = -0.932812808279; 	 
  x[5] = -0.902098806969; 	 
  x[6] = -0.865959503212; 	 
  x[7] = -0.824612230833; 	 
  x[8] = -0.778305651427; 	 
  x[9] = -0.727318255190; 	 
  x[10] = -0.671956684614;  
  x[11] = -0.612553889668;  
  x[12] = -0.549467125095;  
  x[13] = -0.483075801686;  
  x[14] = -0.413779204372;  
  x[15] = -0.341994090826;  
  x[16] = -0.268152185007;  
  x[17] = -0.192697580701;  
  x[18] = -0.116084070675;  
  x[19] = -0.038772417506;  
  x[20] = 0.038772417506;	 
  x[21] = 0.116084070675; 	 
  x[22] = 0.192697580701; 	 
  x[23] = 0.268152185007;	 
  x[24] = 0.341994090826; 	 
  x[25] = 0.413779204372; 	 
  x[26] = 0.483075801686; 	 
  x[27] = 0.549467125095; 	 
  x[28] = 0.612553889668; 	 
  x[29] = 0.671956684614; 	 
  x[30] = 0.727318255190; 	 
  x[31] = 0.778305651427; 	 
  x[32] = 0.824612230833; 	 
  x[33] = 0.865959503212; 	 
  x[34] = 0.902098806969; 	 
  x[35] = 0.932812808279; 	 
  x[36] = 0.957916819214; 	 
  x[37] = 0.977259949984; 	 
  x[38] = 0.990726238699; 	 
  x[39] = 0.998237709711; 	 

  w[0] = 0.004521277099;
  w[1] = 0.010498284521;
  w[2] = 0.016421058382;
  w[3] = 0.022245849194;
  w[4] = 0.027937006980;
  w[5] = 0.033460195283;
  w[6] = 0.038782167974;
  w[7] = 0.043870908186;
  w[8] = 0.048695807635;
  w[9] = 0.053227846984;
  w[10] = 0.057439769099;
  w[11] = 0.061306242493;
  w[12] = 0.064804013457;
  w[13] = 0.067912045815;
  w[14] = 0.070611647391;
  w[15] = 0.072886582396;
  w[16] = 0.074723169058;
  w[17] = 0.076110361901;
  w[18] = 0.077039818164;
  w[19] = 0.077505947978;
  w[20] = 0.077505947978;
  w[21] = 0.077039818164;
  w[22] = 0.076110361901;
  w[23] = 0.074723169058;
  w[24] = 0.072886582396;
  w[25] = 0.070611647391;
  w[26] = 0.067912045815;
  w[27] = 0.064804013457;
  w[28] = 0.061306242493;
  w[29] = 0.057439769099;
  w[30] = 0.053227846984;
  w[31] = 0.048695807635;
  w[32] = 0.043870908186;
  w[33] = 0.038782167974;
  w[34] = 0.033460195283;
  w[35] = 0.027937006980;
  w[36] = 0.022245849194;
  w[37] = 0.016421058382;
  w[38] = 0.010498284521;
  w[39] = 0.004521277099;

  //double sum = 0.;

  for(int i = 0 ; i < 40 ; i++){
    w[i] = (b-a)/2.*w[i];
    x[i] = (b-a)/2.*x[i] + (b+a)/2.;
    //sum += w1*func->Eval(x1);
  }

  //return sum;

}

void GetCoefsIntegrateLaguerre32(std::vector<double> &x, std::vector<double> &w)
{
  // Integrate function from 0 to infinity using Gauss-Laguerre integration
  // with 32 points
  //

  x.resize(32);
  w.resize(32);

  x[0] = 0.044489365833; 	  	
  x[1] = 0.234526109520;	 	
  x[2] = 0.576884629302; 		
  x[3] = 1.072448753818; 		
  x[4] = 1.722408776445; 		
  x[5] = 2.528336706426; 		
  x[6] = 3.492213273022; 		
  x[7] = 4.616456769750; 		
  x[8] = 5.903958504174; 		
  x[9] = 7.358126733186; 		
  x[10] = 8.982940924213; 		
  x[11] = 10.783018632540; 	
  x[12] = 12.763697986743; 	
  x[13] = 14.931139755523; 	
  x[14] = 17.292454336715; 	
  x[15] = 19.855860940336; 	
  x[16] = 22.630889013197; 	
  x[17] = 25.628636022459; 	
  x[18] = 28.862101816323; 	
  x[19] = 32.346629153965; 	
  x[20] = 36.100494805752; 	
  x[21] = 40.145719771539; 	
  x[22] = 44.509207995755; 	
  x[23] = 49.224394987309; 	
  x[24] = 54.333721333397; 	
  x[25] = 59.892509162134; 	
  x[26] = 65.975377287935; 	
  x[27] = 72.687628090663;	
  x[28] = 80.187446977914;	
  x[29] = 88.735340417892; 	
  x[30] = 98.829542868284;	
  x[31] = 111.751398097938; 	

  w[0] = 0.114187105768;
  w[1] = 0.266065216898;
  w[2] = 0.418793137325;
  w[3] = 0.572532846500;
  w[4] = 0.727648788381;
  w[5] = 0.884536719340;
  w[6] = 1.043618875892;
  w[7] = 1.205349274152;
  w[8] = 1.370221338522;
  w[9] = 1.538777256469;
  w[10] = 1.711619352686;
  w[11] = 1.889424063449;
  w[12] = 2.072959340247;
  w[13] = 2.263106633997;
  w[14] = 2.460889072488;
  w[15] = 2.667508126397;
  w[16] = 2.884392092922;
  w[17] = 3.113261327040;
  w[18] = 3.356217692596;
  w[19] = 3.615869856484;
  w[20] = 3.895513044949;
  w[21] = 4.199394104712;
  w[22] = 4.533114978534;
  w[23] = 4.904270287611;
  w[24] = 5.323500972024;
  w[25] = 5.806333214234;
  w[26] = 6.376614674160;
  w[27] = 7.073526580707;
  w[28] = 7.967693509296;
  w[29] = 9.205040331278;
  w[30] = 11.163013090768;
  w[31] = 15.390180415261;

  /*double sum = 0.;

  for(int i = 0 ; i < 32 ; i++){
    sum += w[i]*func->Eval(x[i]);
  }

  return sum;
  */
}

#endif