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

//
//
// BCC paracrystal, powder average
//
// VERY slow, since the function is so ill-behaved and needs LOTS of quadrature
// points. Adaptive methods were even slower and troublesom to converge,
// although in theory they should be a better choice than blindly increasing the number of points.
//
// 150 points seems to give reasonable reproduction of the peak heights in the paper.
// peak locations are correct
// 76 points of quadrature for the smearing is only a guess, it's not been tested yet.
//
// Original implementation - Danilo Pozzo
//		modified and modernized for more efficient integration SRK Nov 2008
//
//REFERENCE 
//Hideki Matsuoka etal. Physical Review B, Vol 36 Num 3, p1754 1987   ORIGINAL PAPER
//Hideki Matsuoka etal. Physical Review B, Vol 41 Num 6, p3854 1990   CORRECTIONS TO PAPER
//
////////////////////////////////////////////////////



Proc PlotBCC_ParaCrystal(num,qmin,qmax)
	Variable num=100, qmin=0.001, qmax=0.7
	Prompt num "Enter number of data points for model: "
	Prompt qmin "Enter minimum q-value (�^-1) for model: " 
	Prompt qmax "Enter maximum q-value (�^-1) for model: "
//
	Make/O/D/n=(num) xwave_BCC_ParaCrystal, ywave_BCC_ParaCrystal
	xwave_BCC_ParaCrystal =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
	Make/O/D coef_BCC_ParaCrystal = {1,220,0.06,40,3e-6,6.3e-6,0.0}
	make/o/t parameters_BCC_ParaCrystal = {"scale","Nearest Neighbor (A)","distortion, g","Sphere Radius (A)","SLD Sphere (A-2)","SLD Solvent (A-2)", "Background (cm-1)"}	
	Edit parameters_BCC_ParaCrystal, coef_BCC_ParaCrystal
	
	Variable/G root:gNordBCC=150
	
	Variable/G root:g_BCC_ParaCrystal
	g_BCC_ParaCrystal := BCC_ParaCrystal(coef_BCC_ParaCrystal, ywave_BCC_ParaCrystal, xwave_BCC_ParaCrystal)
	Display ywave_BCC_ParaCrystal vs xwave_BCC_ParaCrystal
	ModifyGraph marker=29, msize=2, mode=4
	ModifyGraph grid=1,mirror=2
	ModifyGraph log=0
	Label bottom "q (�\\S-1\\M) "
	Label left "I(q) (cm\\S-1\\M)"
	AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
	
	AddModelToStrings("BCC_ParaCrystal","coef_BCC_ParaCrystal","parameters_BCC_ParaCrystal","BCC_ParaCrystal")
//
End

//
//this macro sets up all the necessary parameters and waves that are
//needed to calculate the  smeared model function.
//
//no input parameters are necessary, it MUST use the experimental q-values
// from the experimental data read in from an AVE/QSIG data file
////////////////////////////////////////////////////
// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
Proc PlotSmearedBCC_ParaCrystal(str)								
	String str
	Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
	
	// if any of the resolution waves are missing => abort
	if(ResolutionWavesMissingDF(str))		//updated to NOT use global strings (in GaussUtils)
		Abort
	endif
	
	SetDataFolder $("root:"+str)
	
	// Setup parameter table for model function
	Make/O/D smear_coef_BCC_ParaCrystal = {1,220,0.06,40,3e-6,6.3e-6,0.0}
	make/o/t smear_param_BCC_ParaCrystal = {"scale","Nearest Neighbor (A)","distortion, g","Sphere Radius (A)","SLD Sphere (A-2)","SLD Solvent (A-2)", "Background (cm-1)"}
	Edit smear_param_BCC_ParaCrystal,smear_coef_BCC_ParaCrystal					//display parameters in a table
	
	// output smeared intensity wave, dimensions are identical to experimental QSIG values
	// make extra copy of experimental q-values for easy plotting
	Duplicate/O $(str+"_q") smeared_BCC_ParaCrystal,smeared_qvals
	SetScale d,0,0,"1/cm",smeared_BCC_ParaCrystal
	
	Variable/G gNordBCC = 150	
	Variable/G gs_BCC_ParaCrystal=0
	gs_BCC_ParaCrystal := fSmearedBCC_ParaCrystal(smear_coef_BCC_ParaCrystal,smeared_BCC_ParaCrystal,smeared_qvals)	//this wrapper fills the STRUCT
	
	Display smeared_BCC_ParaCrystal vs smeared_qvals
	ModifyGraph marker=29,msize=2,mode=4
	ModifyGraph log=0
	Label bottom "q (�\\S-1\\M)"
	Label left "I(q) (cm\\S-1\\M)"
	AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
	
	SetDataFolder root:
	AddModelToStrings("SmearedBCC_ParaCrystal","smear_coef_BCC_ParaCrystal","smear_param_BCC_ParaCrystal","BCC_ParaCrystal")
End

// Threaded version
// Threaded XOP = 2.4 s
// non-threaded, non-XOP = 46.8 s
// = x 19.5 speedup !
//
Function BCC_ParaCrystal(cw,yw,xw) : FitFunc
	Wave cw,yw,xw

/////// NO threading /////////
//#if exists("BCC_ParaCrystalX")
//	yw = BCC_ParaCrystalX(cw,xw)
//#else
//	yw = fBCC_ParaCrystal(cw,xw)
//#endif


/// THREADING ///////

//	Variable t1=StopMSTimer(-2)
		
#if exists("BCC_ParaCrystalX")

	Variable npt=numpnts(yw)
	Variable i,nthreads= ThreadProcessorCount
	variable mt= ThreadGroupCreate(nthreads)

	for(i=0;i<nthreads;i+=1)
	//	Print (i*npt/nthreads),((i+1)*npt/nthreads-1)
		ThreadStart mt,i,BCC_ParaCrystal_T(cw,yw,xw,(i*npt/nthreads),((i+1)*npt/nthreads-1))
	endfor

	do
		variable tgs= ThreadGroupWait(mt,100)
	while( tgs != 0 )

	variable dummy= ThreadGroupRelease(mt)
	
#else
		yw = fBCC_ParaCrystal(cw,xw)		//the Igor, non-XOP, non-threaded calculation, messy to make ThreadSafe
#endif

//	Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6



	return(0)
End



// nothing to change here
//
//AAO version, uses XOP if available
// simply calls the original single point calculation with
// a wave assignment (this will behave nicely if given point ranges)
//
// Threaded Version
ThreadSafe Function BCC_ParaCrystal_T(cw,yw,xw,p1,p2) : FitFunc
	Wave cw,yw,xw
	Variable p1,p2

//	Variable t1=StopMSTimer(-2)

#if exists("BCC_ParaCrystalX")
	yw[p1,p2] = BCC_ParaCrystalX(cw,xw)
#else
	yw[p1,p2] = fBCC_ParaCrystal(cw,xw)
#endif

//	Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6

	return(0)
End


//
// unsmeared model calculation
//
Function fBCC_ParaCrystal(w,x) : FitFunc
	Wave w
	Variable x
	
//	 Input (fitting) variables are not used
//	you would give them nice names
	Variable integral,loLim,upLim
	loLim = 0
	upLim = 2*Pi
	
	Variable/G root:gDumY=0		//root:gDumX=0
	

	Variable scale,Dnn,gg,Rad,contrast,background,yy,latticeScale,sld,sldSolv
	scale = w[0]
	Dnn = w[1] //Nearest neighbor distance A
	gg = w[2] //Paracrystal distortion factor
	Rad = w[3] //Sphere radius
	sld = w[4]
	sldSolv = w[5]
	background = w[6] 
	
	contrast = sld - sldSolv
	
	//Volume fraction calculated from lattice symmetry and sphere radius
	latticeScale = 2*(4/3)*pi*(Rad^3)/((Dnn/((3/4)^0.5))^3)

	NVAR/Z nord=root:gNordBCC
	if(NVAR_Exists(nord)!=1)
		nord=20
	endif
	

	integral = IntegrateFn_N(Integrand_BCC_Outer,loLim,upLim,w,x,nord)
	
	integral *= SphereForm_BCC(Rad,contrast,x)*scale*latticeScale
//	integral *= scale			//testing, returns only Z(q)

	integral += background	
	
	Return (integral)
	
End

// the outer integral is also an integral
Function Integrand_BCC_Outer(w,x,dum)
	Wave w
	Variable x,dum
		
	NVAR yy = root:gDumY		
	yy = dum					// save the current dummy yy for use in the inner loop
	Variable retVal,loLim,upLim
	//
	loLim = 0
	upLim = Pi

	NVAR/Z nord=root:gNordBCC
	if(NVAR_Exists(nord)!=1)
		nord=20
	endif

	retVal = IntegrateFn_N(Integrand_BCC_Inner,loLim,upLim,w,x,nord)
	
	return(retVal)
End

//returns the value of the integrand of the inner integral
Function Integrand_BCC_Inner(w,qq,dum)
	Wave w
	Variable qq,dum
	
	NVAR yy = root:gDumY		//use the yy value from the outer loop
	Variable xx,retVal
	xx = dum

	retVal = BCC_Integrand(w,qq,xx,yy)
	
	return(retVal)
End

Function BCC_Integrand(w,qq,xx,yy)
	Wave w
	Variable qq,xx,yy
	
	Variable retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg
	Dnn = w[1] //Nearest neighbor distance A
	gg = w[2] //Paracrystal distortion factor
//	aa = Dnn*((4/3)^0.5)		//Danilo's version (the paper states |bi| = a)
	aa = Dnn
	Da = gg*aa
	
	temp1 = qq*qq*Da*Da
	temp3 = qq*aa	
	
	retVal = BCCeval(xx,yy,temp1,temp3)
	retVal /=4*Pi
	
	return(retVal)
end

Function BCCeval(Theta,Phi,temp1,temp3)
	Variable Theta,Phi,temp1,temp3
	Variable temp6,temp7,temp8,temp9,temp10
	Variable result
	
	temp6 = sin(Theta)
	temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta)
	temp8 = -1*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta)
	temp9 = -1*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta)
	temp10 = exp((-1/8)*temp1*((temp7^2)+(temp8^2)+(temp9^2)))
	result = ((1-(temp10^2))^3)*temp6/((1-2*temp10*cos(0.5*temp3*(temp7))+(temp10^2))*(1-2*temp10*cos(0.5*temp3*(temp8))+(temp10^2))*(1-2*temp10*cos(0.5*temp3*(temp9))+(temp10^2)))
	
	return (result)
end


Function SphereForm_BCC(radius,delrho,x)					
	Variable radius,delrho,x
	
	// variables are:							
	//[2] radius (�)
	//[3] delrho (�-2)
	//[4] background (cm-1)
	
	// calculates scale * f^2/Vol where f=Vol*3*delrho*(sin(qr)-qrcos(qr))/qr^3
	// and is rescaled to give [=] cm^-1
	
	Variable bes,f,vol,f2
	//
	//handle q==0 separately
	If(x==0)
		f = 4/3*pi*radius^3*delrho*delrho*1e8
		return(f)
	Endif
	
	bes = 3*(sin(x*radius)-x*radius*cos(x*radius))/x^3/radius^3
	vol = 4*pi/3*radius^3
	f = vol*bes*delrho		// [=] �
	// normalize to single particle volume, convert to 1/cm
	f2 = f * f / vol * 1.0e8		// [=] 1/cm
	
	return (f2)	
	
End




///////////////////////////////////////////////////////////////
// smeared model calculation
//
Function SmearedBCC_ParaCrystal(s) : FitFunc
	Struct ResSmearAAOStruct &s

//	the name of your unsmeared model (AAO) is the first argument
	Smear_Model_76(BCC_ParaCrystal,s.coefW,s.xW,s.yW,s.resW)

	return(0)
End


/////////////////////////////////////////////////////////////////
//wrapper to calculate the smeared model as an AAO-Struct
// fills the struct and calls the ususal function with the STRUCT parameter
//
// used only for the dependency, not for fitting
//
Function fSmearedBCC_ParaCrystal(coefW,yW,xW)
	Wave coefW,yW,xW
	
	String str = getWavesDataFolder(yW,0)
	String DF="root:"+str+":"
	
	WAVE resW = $(DF+str+"_res")
	
	STRUCT ResSmearAAOStruct fs
	WAVE fs.coefW = coefW	
	WAVE fs.yW = yW
	WAVE fs.xW = xW
	WAVE fs.resW = resW
	
	Variable err
	err = SmearedBCC_ParaCrystal(fs)
	
	return (0)
End