//| This routine creates Gauss-Hermite abscissa and weights
//|	xGH[0],xGH[1],....,xGH[n*(n+1)/2-1]
//|	wGH[0],wGH[1],...,wGH[n*(n+1)/2-1]
//| For example, 
//|	the1st abscissa & weight of 1-terms Gauss-Hermite sum are indexed by 
//|		xGH[0], wGH[0],
//|	the 1st data of 2-term Gauss-Hermite are xGH[1],wGH[1],
//|	the 2nd data of 2-term Gauss-Hermite are xGH[2],wGH[2],
//|	the 1st data of 3-term Gauss-Hermite are xGH[3],wGH[3],
//|		......, 
//|	the j-th data of n-term Gauss-Hermite are
//|		xGH[n*(n-1)/2+j-1], wGH[n*(n-1)/2+j-1]
//|		......, 
//|	the n-th data of n-term Gauss-Hermite are
//|		xGH[n*(n+1)/2-1], wGH[n*(n+1)/2-1]
//|
function GauHerInit(n0)
	variable n0
	variable n,i,its,j,m;
	variable/d z1,z,pp,p3,p2,p1;
	variable/d EPS,PIM4
	variable MAXIT
	EPS=3.0e-14;
	PIM4=0.7511255444649425;
	MAXIT=256;
	if (n0>198)
		print "Too large number !"
		print "Forced to maximum 198"
		n0=198;
	endif;
	Make /D/O /N=(n0*(n0+1)/2) xGH,wGH;
	
	n=1;
	do
		m=floor((n+1)/2);
		i=1;
		do
			if (i==1)
				z=sqrt(2*n+1)-1.85575*(2*n+1)^(-0.16667);
			else
				if (i==2)
					z-=1.14*n^0.426/z;
				else
					if (i==3)
						z=1.86*z+0.86*xGH[n*(n-1)/2];
					else
						if (i==4)
							z=1.91*z+0.91*xGH[n*(n-1)/2+1];
						else
							z=2*z+xGH[n*(n-1)/2+i-3];
						endif;
					endif;
				endif;
			endif;
			its=1;
			do
				p1=PIM4;
				p2=0.0;
				j=1;
				do
					p3=p2;
					p2=p1;
					p1=z*sqrt(2.0/j)*p2-sqrt((j-1)/j)*p3;
					j+=1;
				while (j<=n);
				pp=sqrt(2*n)*p2;
				z1=z;
				z=z1-p1/pp;
				if (abs(z-z1)<=EPS)
					break;
				endif;
				its+=1;
			while(its<=MAXIT);
			if (its>MAXIT)
				print "too many iterations in GauHerInit"
				return 1;
			endif;
			xGH[n*(n-1)/2+i-1]=-z;
			xGH[n*(n-1)/2+n-i]=z;
			wGH[n*(n-1)/2+i-1]=2/(pp*pp)/sqrt(pi);
			wGH[n*(n-1)/2+n-i]=wGH[n*(n-1)/2+i-1];
			i+=1;
		while (i<=m);
		n+=1;
	while (n<=n0);
	return 0;
end;

function TestGauHer()
	wave xGH,wGH;
	variable n0,i,j;
	variable s;
	n0=floor(sqrt(numpnts(xGH)*2));
	i=1;
	do
		j=0;
		s=0;
		do
			s+=wGH[i*(i-1)/2+j];
			j+=1;
		while (j<i);
		print i,s;
		i+=1;
	while (i<n0);
end;