//------------------------------------------------------------------------------
// Coeffitients for numerical integration in range [-1.,1.]
// Int(f(x),x=-1..1) ~ Sum_(i=0..NGauss) f(x_i)*w_i + Delta
// for details see M.Abramovic, I.Stegan, Special Function Handbook
//------------------------------------------------------------------------------
function GaussQuadrature()
{
	this.Const = new Constants();
	
	this.NGauss = -1;
	this.xGauss = [];
	this.wGauss = [];
	
	this.xGauss64 = [
			-0.999305041735772139457,-0.996340116771955279347,-0.991013371476744320739,
			-0.983336253884625956931,-0.973326827789910963742,-0.961008799652053718919,
			-0.946411374858402816062,-0.929569172131939575821,-0.910522137078502805756,
			-0.889315445995114105853,-0.865999398154092819761,-0.840629296252580362752,
			-0.813265315122797559742,-0.783972358943341407610,-0.752819907260531896612,
			-0.719881850171610826849,-0.685236313054233242564,-0.648965471254657339858,
			-0.611155355172393250249,-0.571895646202634034284,-0.531279464019894545658,
			-0.489403145707052957479,-0.446366017253464087985,-0.402270157963991603696,
			-0.357220158337668115950,-0.311322871990210956158,-0.264687162208767416374,
			-0.217423643740007084150,-0.169644420423992818037,-0.121462819296120554470,
			-0.072993121787799039450,-0.024350292663424432509,
			 0.024350292663424432509, 0.072993121787799039450, 0.121462819296120554470,
			 0.169644420423992818037, 0.217423643740007084150, 0.264687162208767416374,
			 0.311322871990210956158, 0.357220158337668115950, 0.402270157963991603696,
			 0.446366017253464087985, 0.489403145707052957479, 0.531279464019894545658,
			 0.571895646202634034284, 0.611155355172393250249, 0.648965471254657339858,
			 0.685236313054233242564, 0.719881850171610826849, 0.752819907260531896612,
			 0.783972358943341407610, 0.813265315122797559742, 0.840629296252580362752,
			 0.865999398154092819761, 0.889315445995114105853, 0.910522137078502805756,
			 0.929569172131939575821, 0.946411374858402816062, 0.961008799652053718919,
			 0.973326827789910963742, 0.983336253884625956931, 0.991013371476744320739,
			 0.996340116771955279347, 0.999305041735772139457];
	this.wGauss64 = [
			0.001783280721696432947, 0.004147033260562467635, 0.006504457968978362856,
			0.008846759826363947723, 0.011168139460131128819, 0.013463047896718642598,
			0.015726030476024719322, 0.017951715775697343085, 0.020134823153530209372,
			0.022270173808383254159, 0.024352702568710873338, 0.026377469715054658672,
			0.028339672614259483228, 0.030234657072402478868, 0.032057928354851553585,
			0.033805161837151609392, 0.035472213256882383811, 0.037055128540240046040,
			0.038550153178615629129, 0.039953741132720341387, 0.041262563242623528610,
			0.042473515123653589007, 0.043583724529323453377, 0.044590558163756563060,
			0.045491627927418144480, 0.046284796581314417296, 0.046968182816210107325,
			0.047540165714830308662, 0.047999388596458307728, 0.048344762234802957170,
			0.048575467441503426935, 0.048690957009139720383,
			0.048690957009139720383, 0.048575467441503426935, 0.048344762234802957170,
			0.047999388596458307728, 0.047540165714830308662, 0.046968182816210107325,
			0.046284796581314417296, 0.045491627927418144480, 0.044590558163756563060,
			0.043583724529323453377, 0.042473515123653589007, 0.041262563242623528610,
			0.039953741132720341387, 0.038550153178615629129, 0.037055128540240046040,
			0.035472213256882383811, 0.033805161837151609392, 0.032057928354851553585,
			0.030234657072402478868, 0.028339672614259483228, 0.026377469715054658672,
			0.024352702568710873338, 0.022270173808383254159, 0.020134823153530209372,
			0.017951715775697343085, 0.015726030476024719322, 0.013463047896718642598,
			0.011168139460131128819, 0.008846759826363947723, 0.006504457968978362856,
			0.004147033260562467635, 0.001783280721696432947];
	this.xGauss40 = [
			-0.998237709710559200350,-0.990726238699457006453,-0.977259949983774262663,-0.957916819213791655805,
			-0.932812808278676533361,-0.902098806968874296728,-0.865959503212259503821,-0.824612230833311663196,
			-0.778305651426519387695,-0.727318255189927103281,-0.671956684614179548379,-0.612553889667980237953,
			-0.549467125095128202076,-0.483075801686178712909,-0.413779204371605001525,-0.341994090825758473007,
			-0.268152185007253681141,-0.192697580701371099716,-0.116084070675255208483,-0.038772417506050821933,
			0.038772417506050821933, 0.116084070675255208483, 0.192697580701371099716, 0.268152185007253681141,
			0.341994090825758473007, 0.413779204371605001525, 0.483075801686178712909, 0.549467125095128202076,
			0.612553889667980237953, 0.671956684614179548379, 0.727318255189927103281, 0.778305651426519387695,
			0.824612230833311663196, 0.865959503212259503821, 0.902098806968874296728, 0.932812808278676533361,
			0.957916819213791655805, 0.977259949983774262663, 0.990726238699457006453, 0.998237709710559200350];
	this.wGauss40 = [
			0.004521277098533191258, 0.010498284531152813615, 0.016421058381907888713, 0.022245849194166957262,
			0.027937006980023401098, 0.033460195282547847393, 0.038782167974472017640, 0.043870908185673271992,
			0.048695807635072232061, 0.053227846983936824355, 0.057439769099391551367, 0.061306242492928939167,
			0.064804013456601038075, 0.067912045815233903826, 0.070611647391286779695, 0.072886582395804059061,
			0.074723169057968264200, 0.076110361900626242372, 0.077039818164247965588, 0.077505947978424811264,
			0.077505947978424811264, 0.077039818164247965588, 0.076110361900626242372, 0.074723169057968264200,
			0.072886582395804059061, 0.070611647391286779695, 0.067912045815233903826, 0.064804013456601038075,
			0.061306242492928939167, 0.057439769099391551367, 0.053227846983936824355, 0.048695807635072232061,
			0.043870908185673271992, 0.038782167974472017640, 0.033460195282547847393, 0.027937006980023401098,
			0.022245849194166957262, 0.016421058381907888713, 0.010498284531152813615, 0.004521277098533191258];
	this.xGauss20 = [
			-0.993128599185094924786,-0.963971927277913791268,-0.912234428251325905868,
			-0.839116971822218823395,-0.746331906460150792614,-0.636053580726515025453,
			-0.510867001950827098004,-0.373706088715419560673,-0.227785851141645078080,
			-0.076526521133497333755, 0.076526521133497333755, 0.227785851141645078080,
			0.373706088715419560673, 0.510867001950827098004, 0.636053580726515025453,
			0.746331906460150792614, 0.839116971822218823395, 0.912234428251325905868,
			0.963971927277913791268, 0.993128599185094924786];
	this.wGauss20 = [
			0.017614007139152118312, 0.040601429800386941331, 0.062672048334109063570,
			0.083276741576704748725, 0.101930119817240435037, 0.118194531961518417312,
			0.131688638449176626898, 0.142096109318382051329, 0.149172986472603746788,
			0.152753387130725850698, 0.152753387130725850698, 0.149172986472603746788,
			0.142096109318382051329, 0.131688638449176626898, 0.118194531961518417312,
			0.101930119817240435037, 0.083276741576704748725, 0.062672048334109063570,
			0.040601429800386941331, 0.017614007139152118312];
	this.xGauss12 = [
			-0.981560634246719,-0.904117256370475,-0.769902674194305,
			-0.587317954286617,-0.367831498998180,-0.125233408511469,
			0.125233408511469, 0.367831498998180, 0.587317954286617,
			0.769902674194305, 0.904117256370475, 0.981560634246719];
	this.wGauss12 = [
			0.047175336386512, 0.106939325995318, 0.160078328543346,
			0.203167426723066, 0.233492536538355, 0.249147045813403,
			0.249147045813403, 0.233492536538355, 0.203167426723066,
			0.160078328543346, 0.106939325995318, 0.047175336386512];
	this.xGauss8 = [
			-0.960289856497536,-0.796666477413627,-0.525532409916329,-0.183434642495650,
			0.183434642495650, 0.525532409916329, 0.796666477413627, 0.960289856497536];
	this.wGauss8 = [
			0.101228536290376, 0.222381034453374, 0.313706645877887, 0.362683783378362,
			0.362683783378362, 0.313706645877887, 0.222381034453374, 0.101228536290376];

this.doGaussQuadrature = function(N)
{
	NGauss = N;
	xGauss = []; 
	wGauss = [];
	//---------------------------------------------------------------------------
	var
		x,dx, xa,xb,
		p1,p2,p,
		eps = 1.e-12,
		xx = [],
		ww = [],
		n,add,i;
	var plm = new LegendrePolynomial();
	//---------------------------------------------------------------------------
	//  Poisk nulei polinoma P_N(x)
	n=0; dx = 2./N; x = 1*dx/2 + this.Const.small_add; 
	i = N/2; p1 = (i%2==0)?1:-1;
	while(n<N/2)
	{
		p2 = plm.Pl0_cos_th(N,x);
		if(p1*p2<0)// there is a zerro! Searching by half-dividing method
		{
			xa = x-dx; xb = x;
			while(true)
			{
				x = 0.5*(1*xa+1*xb); p = plm.Pl0_cos_th(N,x);
				if(Math.abs(p) < eps)
				{
					xx[n]=x;
					if(n>0) dx = 0.2*(xx[n] - xx[n-1]);
					else    dx = 0.4*xx[n];
					n++;
					x = 1*x + 0.5*dx; p1 = plm.Pl0_cos_th(N,x);
					break;
				}
				//``````````````````````````````````````````````````````````````````
				// oshibka v vybore dx!
				if(Math.abs(xa-xb)<eps && Math.abs(p)>this.Const.small_add)
				{
					dx *= 0.1; x = this.Const.zerro-dx; n = 0; p1 = ((N/2)%2==0)?1:-1;
					break;
				}
				// rasstoyanie mezhdu nulymi P_N slishkom malo, avariinoe zavershenie
				else if(Math.abs(xa-xb)<eps && Math.abs(p)<eps)
				{
					console.log("ERROR in Gauss quadrature search:\n"+
						"\tYou choose too high order of Legendre polynomial!");
					return false;
				}
				//``````````````````````````````````````````````````````````````````
				if(p*p1 < 0) xb=x; else xa=x;
			}
		}
		x += dx;
	}
	//---------------------------------------------------------------------------
	//  Vychislenie coefficientov
	for(n=0;n<N/2;n++)
	{
		p2 = plm.Pl0_cos_th(N,xx[n]);
		p = N*(xx[n]*p2 - plm.PL_1); p = p*p;
		ww[n] = 2.*(1-xx[n]*xx[n])/p;
	}
	//---------------------------------------------------------------------------
	//  Zapolnyaem osnovnye massivy
	if(N%2==0)// if N is odd then P(x=0)=0
	{
		this.xGauss[N/2] = 0; 
		plm.Pl0_cos_th(N,0.0);
		p = plm.PL_1*N; p = p*p; 
		this.wGauss[N/2] = 2./p;
		add = 1;
	} else add=0;
	for(n=0;n<N/2;n++)
	{
		this.xGauss[N/2-1-n] = -xx[n]; this.xGauss[N/2+add+n] = xx[n];
		this.wGauss[N/2-1-n] =  ww[n]; this.wGauss[N/2+add+n] = ww[n];
	}
	//---------------------------------------------------------------------------
	return true;
}
this.setGaussQuadrature = function(N)
{
    switch(N)
    {
        default:
            if(N>0) return this.doGaussQuadrature(N);
            return false;
        case  8:
            this.NGauss=N; this.xGauss = this.xGauss8; this.wGauss = this.wGauss8; return true;
        case 12:
            this.NGauss=N; this.xGauss = this.xGauss12; this.wGauss = this.wGauss12; return true;
        case 20:
            this.NGauss=N; this.xGauss = this.xGauss20; this.wGauss = this.wGauss20; return true;
        case 40:
            this.NGauss=N; this.xGauss = this.xGauss40; this.wGauss = this.wGauss40; return true;
        case 64:
            this.NGauss=N; this.xGauss = this.xGauss64; this.wGauss = this.wGauss64; return true;
    }
}
}