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

////////////////////////////////////////////////////
//
// A template for writing your own scattering functions
//
// FIVE things to do:
// 1) identify the equation and the free variables before coding
//
// --then you may start coding ---
//
// 2) change coef_MyModel to have the coefficients (defined in w[])
//    in the exact order you have specified
//    - this change is one line, search for "CH#1"
// 3) change parameters_MyModel to has parameter names, that match
//    the coefficients - this line follows the coefficients of step 1
//	    - this is "CH#2"
// 4) change the function MyModel() to calculate what you want
//    w[N] must contain all the fitting parameters - search for "CH#3"
// 5) repeat steps 2 and 3 (copy/paste) changing the RHS of the same
//    lines in the macro for the smeared model calculation
//    - search for "CH#4", and the following line
//
// - then, all you need to do is (re)plot the macro from the Macros menu
//
// As-is, the model calculates a form factor for a uniform sphere 
//
// v 2.00 -- 01 MAY 2006 SRK
// v 2.1 -- Jul 2007 SRK
//				converted to use structures, AAO, like the rest of the package
//
////////////////////////////////////////////////////

//this macro sets up all the necessary parameters and waves that are
//needed to calculate the model function.
//
Macro PlotMyModel(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_MyModel, ywave_MyModel
	xwave_MyModel =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
	Make/O/D coef_MyModel = {0.01,60,6e-6,0.01}			//CH#1
	make/o/t parameters_MyModel = {"scale","Radius (A)","contrast (�-2)","bkgd (cm-1)"}	//CH#2
	Edit parameters_MyModel, coef_MyModel
	
	Variable/G root:g_MyModel
	g_MyModel := MyModel(coef_MyModel, ywave_MyModel, xwave_MyModel)
	Display ywave_MyModel vs xwave_MyModel
	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("MyModel","coef_MyModel","MyModel")
//
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
Macro PlotSmearedMyModel(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_MyModel = {0.01,60,6e-6,0.0}		//CH#4 change this line and the following line
	make/o/t smear_parameters_MyModel = {"scale","Radius (A)","contrast (�-2)","bkgd (cm-1)"}
	Edit smear_parameters_MyModel,smear_coef_MyModel					//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_MyModel,smeared_qvals
	SetScale d,0,0,"1/cm",smeared_MyModel
					
	Variable/G gs_MyModel=0
	gs_MyModel := fSmearedMyModel(smear_coef_MyModel,smeared_MyModel,smeared_qvals)	//this wrapper fills the STRUCT
	
	Display smeared_MyModel 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("SmearedMyModel","smear_coef_MyModel","MyModel")
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)
Function MyModel(cw,yw,xw) : FitFunc
	Wave cw,yw,xw
	
#if exists("MyModelX")
	yw = MyModelX(cw,xw)
#else
	yw = fMyModel(cw,xw)
#endif
	return(0)
End


//CH#3
// you should write your function to calculate the intensity
// for a single q-value (that's the input parameter x)
// based on the wave (array) of parameters that you send it (w)
//
// unsmeared model calculation
//
Function fMyModel(w,x) : FitFunc
	Wave w
	Variable x
	
	//	 Input (fitting) variables are:						
	//[0] scale
	//[1] radius (�)
	//[2] delrho (�-2)
	//[3] background (cm-1)
	
	Variable scale,radius,delrho,bkg,qval			
	scale = w[0]
	radius = w[1]
	delrho = w[2]
	bkg = w[3]
	
	qval = x		//rename the input q-value, purely for readability
	
	// 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 to avoid divide-by-zero in the Bessel func
	If(qval==0)
		f = 4/3*pi*radius^3*delrho*delrho*scale*1e8 + bkg
		return(f)
	Endif
	
	bes = 3*(sin(qval*radius)-qval*radius*cos(qval*radius))/qval^3/radius^3
	vol = 4*pi/3*radius^3
	f = vol*bes*delrho		// [=] �
	// normalize to single particle volume and convert to 1/cm
	f2 = f * f / vol * 1e8		// [=] 1/cm
	
	return (scale*f2+bkg)	// Scale, then add in the background
	
End

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

//	the name of your unsmeared model (AAO) is the first argument
	Smear_Model_20(MyModel,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 fSmearedMyModel(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 = SmearedMyModel(fs)
	
	return (0)
End