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

//-------------------------------------------------------------------------
//|
//| Function A_cylinder(muR, degTheta) 
//| Function A_sphere(muR, degTheta)
//|
//| Overall transmittance of cylinder and sphere
//| as functions of muR & theta.
//|
//| The algorithm has been proposed by the author 
//| [T. Ida, J. Appl. Cryst. 43, 1124-1125 (2010)].
//|
//| Absorption correction factors will be given by 
//| "1/A_cylinder(muR, degTheta)" or "1/A_sphere(muR, degTheta)".
//| See Table 6.3.3.2 & 6.3.3.3 (p. 596) in 
//| "International Tables for Crytallography Vol. C".
//| Note that the values calculated by this method are more accurate 
//| than the values in the "International Tables".
//| The mathematical formulas are given in a literature
//| [G. Thorkildsen & H. B. Larsen, Acta Cryst. A54, 186-190 (1998)]. 
//|
//------------------------------------------------------------------------- 
function A_cylinder(muR, degTheta)
	variable muR;
		// product of linear absorption coefficient (mu)
		// and radius of cylinder (R)
	variable degTheta;	// Bragg angle in [deg.]

	variable/g gGauLegInitialized;
	variable theta;	// Bragg angle in [rad.]

	variable nX = 12, nY = 12;
	variable iX, iY;
	variable x0, x1, x;
	variable y0, y1, y;
	variable weightX, weightY;
	variable a, h;
	variable ans;

	make/d/o/n=1 wA={2};	// parameter used in HyperGPFQ()
	make/d/o/n=2 wB={1.5, 2.5};	// parameters used in HyperGPFQ()

	if (muR == 0)
		return 1.0;
	endif;
	if (degTheta == 0)
		ans = 2 * BessI(2, 2*muR);
		ans += BessI(1, 2*muR) / muR;
		ans -= 16*muR/(3*pi) * HyperGPFQ(wA, wB, muR*muR);
		return ans;
	endif;
	if (degTheta == 90)
		ans = BessI(1, 4*muR) / (2*muR);
		ans -= 16*muR/(3*pi)*HyperGPFQ(wA, wB, 4*muR*muR);
		return ans;
	endif;

	if (gGauLegInitialized	== 0)
		GauLegInit(128);
		gGauLegInitialized = 1;
	endif;
	wave xGL, wGL;

	theta = degTheta * pi / 180;
	x0 = 0;
	x1 = pi - 2*theta;
	y0 = 0;
	y1 = 2*theta;
	ans = 0;
	iX = 0;
	do
		x = x0 + (x1 - x0) * xGL[nX*(nX-1)/2+iX];
		weightX = wGL[nX*(nX-1)/2+iX];
		iY = 0;
		do
			y = y0 + (y1 - y0) * xGL[nY*(nY-1)/2+iY];
			weightY = wGL[nY*(nY-1)/2+iY];
			a = 2 * muR * sin(x) * cos(y - theta) / cos(theta);
			h = sin(x+y) * sin(x-y+2*theta) * exp(-a);
			ans += weightX * weightY * h;
			iY += 1;
		while (iY < nY);
		iX += 1;
	while (iX < nX);
	ans *= (x1 - x0) * (y1 - y0);
	ans *= 2 / (pi*sin(2*theta));
	return ans;
end function;

function A_sphere(muR, degTheta)
	variable muR;	// product of linear absorption coefficient and radius of sphere
	variable degTheta;	// Bragg angle in [deg.]

	variable/g gGauLegInitialized;
	variable theta;	// Bragg angle in [rad.]

	variable nX = 12, nY = 12;
	variable iX, iY;
	variable x0, x1, x;
	variable y0, y1, y;
	variable weightX, weightY;
	variable a, h;
	variable ans;

	if (muR == 0)
		return 1.0;
	endif;
	if (degTheta == 0)
		ans = 1 - (1 + 2*muR + 2*muR^2) * exp(-2*muR);
		ans *= 3/(4*muR^3);
		return ans;
	endif;
	if (degTheta == 90)
		ans = 1 - exp(4*muR) + 4*muR + 8*exp(4*muR)*muR^2;
		ans *= 3/(64*muR^3) * exp(-4*muR);
		return ans;
	endif;

	if (gGauLegInitialized	== 0)
		GauLegInit(128);
		gGauLegInitialized = 1;
	endif;
	wave xGL, wGL;
	make/d/o/n=1 wA={2};	// parameter used in HyperGPFQ()
	make/d/o/n=2 wB={0.5, 2.5};	// parameters used in HyperGPFQ()

	theta = degTheta * pi / 180;
	x0 = 0;
	x1 = pi - 2*theta;
	y0 = 0;
	y1 = 2*theta;
	ans = 0;
	iX = 0;
	do
		x = x0 + (x1 - x0) * xGL[nX*(nX-1)/2+iX];
		weightX = wGL[nX*(nX-1)/2+iX];
		iY = 0;
		do
			y = y0 + (y1 - y0) * xGL[nY*(nY-1)/2+iY];
			weightY = wGL[nY*(nY-1)/2+iY];
			a = 2 * muR * sin(x) * cos(y - theta) / cos(theta);
			h = HyperGPFQ(wA, wB, a*a/4);
			h -= 9*pi/(4*a)*BesselI(2, a);
			h -= 3*pi/4 * BesselI(3, a);
			h *= sin(x+y) * sin(x-y+2*theta);
			ans += weightX * weightY * h;
			iY += 1;
		while (iY < nY);
		iX += 1;
	while (iX < nX);
	ans *= (x1 - x0) * (y1 - y0);
	ans *= 2 / (pi*sin(2*theta));
	return ans;
end function;

//--------------------------------------------------------------------
//| The following 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 1st abscissa & weight of 1-term Gauss-Legendre sum is 
//| 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]
//--------------------------------------------------------------------
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;