//| This routine creates Gauss-Legendre abscissa and weights
//| for 0<x<1 region as
//|	xGL[0],xGL[1],....,xGL[n*(n+1)/2-1]
//|	wGL[0],wGL[1],...,wGL[n*(n+1)/2-1]
//| For example, 
//|	the1st abscissa & weight of 1-terms Gauss-Legendre sum are indexed by 
//|		xGL[0], wGL[0],
//|	the 1st data of 2-term Gauss-Legendre are xGL[1],wGL[1],
//|	the 2nd data of 2-term Gauss-Legendre are xGL[2],wGL[2],
//|	the 1st data of 3-term Gauss-Legendre are xGL[3],wGL[3],
//|		......, 
//|	the j-th data of n-term Gauss-Legendre are
//|		xGL[n*(n-1)/2+j-1], wGL[n*(n-1)/2+j-1]
//|		......, 
//|	the n-th data of n-term Gauss-Legendre are
//|		xGL[n*(n+1)/2-1], wGL[n*(n+1)/2-1]
//|
function GauLegInit(n0)
	variable n0
	variable/d x1,x2;
	variable m,n,nL,j,i;
	variable/d z1,z,xm,xl,pp,p3,p2,p1;

	x1=0;
	x2=1;
	Make /D/O /N=(n0*(n0+1)/2) xGL,wGL;
	
	xGL[0]=0.5;
	wGL[0]=1;
	nL=1;
	n=2;
	do
		m=floor((n+1)/2);
		xm=0.5*(x2+x1);
		xl=0.5*(x2-x1);
		i=1;
		do
			z=cos(Pi*(i-0.25)/(n+0.5));
			do
				p1=1.0;
				p2=0.0;
				j=1;
				do
					p3=p2;
					p2=p1;
					p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
					j+=1;
				while (j<n+0.5);
				pp=n*(z*p1-p2)/(z*z-1.0);
				z1=z;
				z=z1-p1/pp;
			while (abs(z-z1) > 3.0e-11);
			xGL[nL+i-1]=xm-xl*z;
			xGL[nL+n-i]=xm+xl*z;
			wGL[nL+i-1]=2.0*xl/((1.0-z*z)*pp*pp);
			wGL[nL+n-i]=wGL[nL+i-1];
			i+=1;
		while (i<m+0.5);
		nL+=n;
		n+=1;
	while (n<=n0);
	return 0;
end;