#pragma rtGlobals=1		// Use modern global access method.

//---------------------------------------------------------------------
//| Transmittance of cylinder, 
//| based on the formulae proposed by Thorkildsen & Larsen
//| [ G. Thorkildsen & H. B. Larsen, Acta Cryst. A54, 172-185 (1998);
//|   G. Thorkildsen & H. B. Larsen, Acta Cryst. A54, 186-190 (1998). ]
//| Efficiency of combination with Gauss-Legendre quadrature 
//| has been confirmed by Ida
//| [ T. Ida, J. Appl. Cryst. 43, 1124-1125 (2010) ]
//|
//| uses static function fThorkildsen(muR, theta, x, y)
//| uses static function AThorkildsen90(muR, nY0)
//| uses static function GauLegInit(n0)
//|
//| To test the values of A* = 1/A (absorption correction factor)
//| [ Internatinal Tables of Crystallography Vol. C p. 596 ] 
//| for muR = 0.5, theta = 30 deg., for example, 
//| type "print 1/AThorkildsen(0.5, 30*pi/180, 8, 8);"
//| on command-line.
//---------------------------------------------------------------------
function AThorkildsen(muR, theta, nX0, nY0)
	variable muR;	// mu * R, 
		// mu : linear absorption coefficient
		// R : radius of cylinder
	variable theta;	// Bragg angle in [rad.]
	variable nX0;	// number of sampling points about x, typically = 8
	variable nY0;	// number of sampling points about y, typically = 8

	nvar gAlgorithm;
	gAlgorithm = 2;	// method of numerical integral
		// = 0: Mid-point, 1: Simpson, 2: Gauss-Legendre

	if ((exists("xGL") == 0) %| (exists("wGL") == 0))
		GauLegInit(128);	// -> "GaussLegendreProc"; 
			// initialize waves xGL, wGL for Gauss-Legendre method
	endif;
	wave xGL, wGL;	// Gauss-Legendre abscissa & weights

	variable a, ans;
	variable alphaX, betaX, x;
	variable alphaY, betaY, y;
	variable iX, iY, nX, nY;
	variable wX, wY;

	if (theta==pi/2)
		return AThorkildsen90(muR,nY0);
	endif;
	if (gAlgorithm == 1)
		nX = 2*ceil(nX0/2);
		nY = 2*ceil(nY0/2);
	else
		nX = nX0;
		nY = nY0;
	endif;

	ans = 0;
	alphaY = 0;
	betaY = 2*theta;
	iY = 0;
	do
		if (gAlgorithm == 1)	// Simspson's method
			y = alphaY + (betaY-alphaY)*iY/nY;
			wY = ((mod(iY,nY)==0) ? (1/3) : ((mod(iY,2)==1) ? (4/3) : (2/3)));
			wY = wY / nY;
		elseif (gAlgorithm == 2)	// Gauss-Legendre method
			y = alphaY + (betaY-alphaY) * xGL[nY*(nY-1)/2+iY];
			wY = ((iY==nY) ? (0) : (wGL[nY*(nY-1)/2+iY]));
		else	// mid-point method
			y = alphaY + (betaY-alphaY)*(iY+0.5)/nY;
			wY = ((iY==nY) ? (0) : (1/nY));
		endif;
		alphaX = 0;
		betaX = pi - 2*theta;
		iX = 0;
		do
			if (gAlgorithm == 1)	// Simspson's method
				x = alphaX + (betaX-alphaX)*iX/nX;
				wX = ((mod(iX,nX)==0) ? (1/3) : ((mod(iX,2)==1) ? (4/3) : (2/3)));
				wX = wX/nX;
			elseif (gAlgorithm == 2)	// Gauss-Legendre method
				x = alphaX + (betaX-alphaX)*xGL[nX*(nX-1)/2+iX];
				wX = ((iX==nX) ? (0) : (wGL[nX*(nX-1)/2+iX]));
			else	// mid-point method
				x = alphaX + (betaX-alphaX)*(iX+0.5)/nX;
				wX = ((iX==nX) ? (0) : (1/nX));
			endif;
			
			ans += wX * wY * fThorkildsen(muR,theta,x,y);

			iX += 1;
		while(iX < nX+1);
		iY += 1;
	while(iY < nY+1);
	if (theta==0)
		ans *= 2/pi*(betaX-alphaX);
	else
		ans *= 2/(pi*sin(2*theta))*(betaX-alphaX)*(betaY-alphaY);
	endif;
	return ans;
end;

//----------------------------------------------
//| Integrand function, called by AThorkildsen()
//----------------------------------------------
static function fThorkildsen(muR, theta, x, y)
	variable muR, theta, x, y;
	variable a, ans;
	a = 2 * muR * sin(x) * cos(y-theta) / cos(theta);
	ans = sin(x+y) * sin(x-y+2*theta) * exp(-a);
	return ans;
end;

//------------------------------------------------------------
//| Special case for Theta = 90 deg., called by AThorkildsen()
//------------------------------------------------------------
static function AThorkildsen90(muR, nY0)
	variable muR, nY0;
	nvar gAlgorithm;	// = 0: Mid-point, 1: Simpson, 2: Gauss-Legendre
	wave xGL, wGL;
	variable a, ans;
	variable alphaY, betaY, y, wY;
	variable iY, nY;
	if (gAlgorithm == 1)
		nY = 2*ceil(nY0/2);
	else
		nY = nY0;
	endif;
	ans = 0;
	alphaY = 0;
	betaY = pi;
	iY = 0;
	do
		if (gAlgorithm == 1)	// Simspson's method
			y = alphaY + (betaY-alphaY)*iY/nY;
			wY = ((mod(iY,nY)==0) ? (1/3) : ((mod(iY,2)==1) ? (4/3) : (2/3)));
			wY=wY/nY;
		elseif (gAlgorithm == 2)	// Gauss-Legendre method
			y = alphaY + (betaY-alphaY)*xGL[nY*(nY-1)/2+iY];
			wY = ((iY == nY) ? (0) : (wGL[nY*(nY-1)/2+iY]));
		else	// mid-point method
			y = alphaY + (betaY-alphaY)*(iY+0.5)/nY;
			wY = ((iY==nY) ? (0) : (1/nY));
		endif;
		ans += wY*(1-exp(-4*muR*sin(y)))*sin(y);
		iY += 1;
	while(iY < nY+1);
	ans *= (betaY-alphaY) / (2*pi*muR);
	return ans;
end;
			
//---------------------------------------------------------------------
//| 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, 
//|	the 0th-abscissa & weight of 1-term Gauss-Legendre sum 
//|		are indexed by xGL[0], wGL[0],
//|
//|	the 0th data of 2-term Gauss-Legendre are xGL[1],wGL[1],
//|	the 1st data of 2-term Gauss-Legendre are xGL[2],wGL[2],
//|
//|	the 0th data of 3-term Gauss-Legendre are xGL[3],wGL[3],
//|
//|		......, 
//|
//|	the j-th (zero-base) data of n-term Gauss-Legendre are
//|		xGL[n*(n-1)/2+j], wGL[n*(n-1)/2+j]
//|
//|		......, 
//|
//|	the (n-1)-th (zero-base) data of n-term Gauss-Legendre are
//|		xGL[n*(n+1)/2-1], wGL[n*(n+1)/2-1]
//---------------------------------------------------------------------
static 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;