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

////////////////////////////////////////////////////
//
// paracrystalline model of Lamellar stacks from J. Pedersen
//
// non-integer # of stacks are calculated as linear combination of trunc(N) and trunc(N+1) 
//
// References:
// - J. Appl. cryst. 30 (1997) 975.
// - JPC B 103 (1999) 9888-9897.
//
////////////////////////////////////////////////////

//
Proc PlotLamellar_ParaCrystal(num,qmin,qmax)
	Variable num=200, 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_LamParaCryst, ywave_LamParaCryst
	xwave_LamParaCryst =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
	Make/O/D coef_LamParaCryst = {1,33,20,250,0.2,1e-6,6.34e-6,0}		
	make/o/t parameters_LamParaCryst = {"scale","Lamellar thickness (A)","N Layers","layer spacing (A)","polydisp of spacing","SLD of Layer (A-2)","SLD of solvent (A-2)","Incoherent Bgd (cm-1)"}	//CH#2
	Edit parameters_LamParaCryst, coef_LamParaCryst
	
	Variable/G root:g_LamParaCryst
	g_LamParaCryst := Lamellar_ParaCrystal(coef_LamParaCryst, ywave_LamParaCryst, xwave_LamParaCryst)
	Display ywave_LamParaCryst vs xwave_LamParaCryst
	ModifyGraph marker=29, msize=2, mode=4
	ModifyGraph log=1,grid=1,mirror=2
	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("Lamellar_ParaCrystal","coef_LamParaCryst","parameters_LamParaCryst","LamParaCryst")
//
End


//
//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 PlotSmearedLamellar_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_LamParaCryst = {1,33,20,250,0.2,1e-6,6.34e-6,0}		
	make/o/t smear_parameters_LamParaCryst = {"scale","Lamellar thickness (A)","N Layers","layer spacing (A)","polydisp of spacing","SLD of Layer (A-2)","SLD of solvent (A-2)","Incoherent Bgd (cm-1)"}
	Edit smear_parameters_LamParaCryst,smear_coef_LamParaCryst					//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_LamParaCryst,smeared_qvals
	SetScale d,0,0,"1/cm",smeared_LamParaCryst
					
	Variable/G gs_LamParaCryst=0
	gs_LamParaCryst := fSmearedLamellar_ParaCrystal(smear_coef_LamParaCryst,smeared_LamParaCryst,smeared_qvals)	//this wrapper fills the STRUCT
	
	Display smeared_LamParaCryst vs smeared_qvals
	ModifyGraph log=1,marker=29,msize=2,mode=4
	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("SmearedLamellar_ParaCrystal","smear_coef_LamParaCryst","smear_parameters_LamParaCryst","LamParaCryst")
End


//
//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)
Function Lamellar_ParaCrystal(cw,yw,xw) : FitFunc
	Wave cw,yw,xw
	
#if exists("Lamellar_ParaCrystalX")
	yw = Lamellar_ParaCrystalX(cw,xw)
#else
	yw = fLamellar_ParaCrystal(cw,xw)
#endif
	return(0)
End

//
// unsmeared model calculation
//
Function fLamellar_ParaCrystal(w,x) : FitFunc
	Wave w
	Variable x
	
//	 Input (fitting) variables are:
	//[0] scale factor
	//[1]	thickness
	//[2]	number of layers
	//[3]	spacing between layers
	//[4]	polydispersity of spacing
	//[5] SLD lamellar
	//[6] SLD solvent
	//[7] incoherent background
//	give them nice names
	Variable inten,qval,scale,th,nl,davg,pd,contr,bkg
	Variable xi,ww,Pbil,Zq,Znq,Snq,an,sldLayer,sldSolvent
	
	scale = w[0]
	th = w[1]
	nl = w[2]
	davg = w[3]
	pd = w[4]
	sldLayer = w[5]
	sldSolvent = w[6]
	bkg = w[7]
	
	contr = w[5] - w[6]
	qval = x
		
	//get the fractional part of nl, to determine the "mixing" of N's
	Variable n1,n2,xn
	
	n1 = trunc(nl)		//rounds towards zero
	n2 = n1 + 1
	xn = n2 - nl			//fractional contribution of n1
	
	ww = exp(-qval*qval*pd*pd*davg*davg/2)
	
	//calculate the n1 contribution
	an = paraCryst_an(ww,qval,davg,n1)
	Snq = paraCryst_sn(ww,qval,davg,n1,an)
	
	Znq = xn*Snq
	
	//calculate the n2 contribution
	an = paraCryst_an(ww,qval,davg,n2)
	Snq = paraCryst_sn(ww,qval,davg,n2,an)
	
	Znq += (1-xn)*Snq
	
	//and the independent contribution
	Znq += (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg))
	
	//the limit when NL approaches infinity
//	Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg))
	
	xi = th/2		//use 1/2 the bilayer thickness
	Pbil = (sin(qval*xi)/(qval*xi))^2
	
	inten = 2*pi*contr*contr*Pbil*Znq/qval^2
	inten *= 1e8
	
	Return (scale*inten+bkg)
End

Function paraCryst_sn(ww,qval,davg,nl,an)
	Variable ww,qval,davg,nl,an
	
	Variable snq
	Snq = an/(nl*(1+ww^2-2*ww*cos(qval*davg))^2)
	return(snq)
end


Function paraCryst_an(ww,qval,davg,nl)
	Variable ww,qval,davg,nl
	
	Variable an
	an = 4*ww^2 - 2*(ww^3+ww)*cos(qval*davg)
	an -= 4*ww^(nl+2)*cos(nl*qval*davg)
	an += 2*ww^(nl+3)*cos((nl-1)*qval*davg)
	an += 2*ww^(nl+1)*cos((nl+1)*qval*davg)
	return(an)
end


///////////////////////////////////////////////////////////////
// smeared model calculation
//
// you don't need to do anything with this function, as long as
// your Lamellar_ParaCrystal works correctly, you get the resolution-smeared
// version for free.
//
// this is all there is to the smeared model calculation!
Function SmearedLamellar_ParaCrystal(s) : FitFunc
	Struct ResSmearAAOStruct &s

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

	return(0)
End


///////////////////////////////////////////////////////////////


// nothing to change here
//
//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 fSmearedLamellar_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 = SmearedLamellar_ParaCrystal(fs)
	
	return (0)
End