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

////////////////////////////////////////////////
//
// this function is for the form factor of a sphere with some
// number of shells around a central core (currently 1-2-3)
//
// monodisperse and polydisperse (and smeared) versions are included
// - for the polydisperse models, only polydispersity of the core is taken
// into account, and done numerically. for a Schulz distribution, this
// should be possible to do analytically, whith a great savings in computation
// time.
//
// It may also be useful to think of scenarios where the layers as well are
// polydisperse - to break up the very regular spacing of the layers, which
// is not a very natural structure.
//
// 03 MAR 04 SRK
////////////////////////////////////////////////
// Four shell model added and polyCore versions written in separate ipf (polyCore_and_NShells.ipf)
//
// Four shell model added at request of H Wacklin to model lipid vesicles with varying deuteration
//
// DEC 2008 AJJ
//////////////////////////////////////////////

//this macro sets up all the necessary parameters and waves that are
//needed to calculate the model function.
//
Proc PlotOneShell(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_OneShell, ywave_OneShell
	xwave_OneShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
	Make/O/D coef_OneShell = {1.,60,6.4e-6,10,1e-6,6.4e-6,0.001}
	make/o/t parameters_OneShell = {"scale","core radius (A)","Core SLD (A-2)","Shell thickness (A)","Shell SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
	Edit parameters_OneShell, coef_OneShell
	
	Variable/G root:g_OneShell
	g_OneShell := OneShell(coef_OneShell, ywave_OneShell, xwave_OneShell)
	Display ywave_OneShell vs xwave_OneShell
	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("OneShell","coef_OneShell","parameters_OneShell","OneShell")
//
End

Proc PlotTwoShell(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_TwoShell, ywave_TwoShell
	xwave_TwoShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
	Make/O/D coef_TwoShell = {1.,60,6.4e-6,10,1e-6,10,2e-6,6.4e-6,0.001}
	make/o/t parameters_TwoShell = {"scale","core radius (A)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
	Edit parameters_TwoShell, coef_TwoShell
	
	Variable/G root:g_TwoShell
	g_TwoShell := TwoShell(coef_TwoShell, ywave_TwoShell, xwave_TwoShell)
	Display ywave_TwoShell vs xwave_TwoShell
	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("TwoShell","coef_TwoShell","parameters_TwoShell","TwoShell")
//
End

Proc PlotThreeShell(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_ThreeShell, ywave_ThreeShell
	xwave_ThreeShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
	Make/O/D coef_ThreeShell ={1.,60,6.4e-6,10,1e-6,10,2e-6,10,3e-6,6.4e-6,0.001}
	make/o/t parameters_ThreeShell = {"scale","core radius (A)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Shell 3 thickness","Shell 3 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
	Edit parameters_ThreeShell, coef_ThreeShell
	
	Variable/G root:g_ThreeShell
	g_ThreeShell := ThreeShell(coef_ThreeShell, ywave_ThreeShell, xwave_ThreeShell)
	Display ywave_ThreeShell vs xwave_ThreeShell
	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("ThreeShell","coef_ThreeShell","parameters_ThreeShell","ThreeShell")
//
End

Proc PlotFourShell(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_FourShell, ywave_FourShell
	xwave_FourShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
	Make/O/D coef_FourShell ={1.,60,6.4e-6,10,1e-6,10,2e-6,10,3e-6,10,4e-6,6.4e-6,0.001}
	make/o/t parameters_FourShell = {"scale","core radius (A)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Shell 3 thickness","Shell 3 SLD (A-2)","Shell 4 thickness","Shell 4 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
	Edit parameters_FourShell, coef_FourShell
	
	Variable/G root:g_FourShell
	g_FourShell := FourShell(coef_FourShell, ywave_FourShell, xwave_FourShell)
	Display ywave_FourShell vs xwave_FourShell
	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("FourShell","coef_FourShell","parameters_FourShell","FourShell")
//
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 PlotSmearedOneShell(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_OneShell =  {1.,60,6.4e-6,10,1e-6,6.4e-6,0.001}
	make/o/t smear_parameters_OneShell =  {"scale","core radius (A)","Core SLD (A-2)","Shell thickness (A)","Shell SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
	Edit smear_parameters_OneShell,smear_coef_OneShell					//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_OneShell,smeared_qvals
	SetScale d,0,0,"1/cm",smeared_OneShell
					
	Variable/G gs_OneShell=0
	gs_OneShell := fSmearedOneShell(smear_coef_OneShell,smeared_OneShell,smeared_qvals)	//this wrapper fills the STRUCT
	
	Display smeared_OneShell 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("SmearedOneShell","smear_coef_OneShell","smear_parameters_OneShell","OneShell")
End

Proc PlotSmearedTwoShell(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_TwoShell =  {1.,60,6.4e-6,10,1e-6,10,2e-6,6.4e-6,0.001}
	make/o/t smear_parameters_TwoShell =  {"scale","core radius (A)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
	Edit smear_parameters_TwoShell,smear_coef_TwoShell					//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_TwoShell,smeared_qvals
	SetScale d,0,0,"1/cm",smeared_TwoShell
					
	Variable/G gs_TwoShell=0
	gs_TwoShell := fSmearedTwoShell(smear_coef_TwoShell,smeared_TwoShell,smeared_qvals)	//this wrapper fills the STRUCT
	
	Display smeared_TwoShell 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("SmearedTwoShell","smear_coef_TwoShell","smear_parameters_TwoShell","TwoShell")
End


Proc PlotSmearedThreeShell(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_ThreeShell =  {1.,60,6.4e-6,10,1e-6,10,2e-6,10,3e-6,6.4e-6,0.001}
	make/o/t smear_parameters_ThreeShell =  {"scale","core radius (A)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Shell 3 thickness","Shell 3 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
	Edit smear_parameters_ThreeShell,smear_coef_ThreeShell					//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_ThreeShell,smeared_qvals
	SetScale d,0,0,"1/cm",smeared_ThreeShell
					
	Variable/G gs_ThreeShell=0
	gs_ThreeShell := fSmearedThreeShell(smear_coef_ThreeShell,smeared_ThreeShell,smeared_qvals)	//this wrapper fills the STRUCT
	
	Display smeared_ThreeShell 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("SmearedThreeShell","smear_coef_ThreeShell","smear_parameters_ThreeShell","ThreeShell")
End

Proc PlotSmearedFourShell(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_FourShell = {1.,60,6.4e-6,10,1e-6,10,2e-6,10,3e-6,10,4e-6,6.4e-6,0.001}
	make/o/t smear_parameters_FourShell =  {"scale","core radius (A)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Shell 3 thickness","Shell 3 SLD (A-2)","Shell 4 thickness","Shell 4 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"}
	Edit smear_parameters_FourShell,smear_coef_FourShell					//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_FourShell,smeared_qvals
	SetScale d,0,0,"1/cm",smeared_FourShell
					
	Variable/G gs_FourShell=0
	gs_FourShell := fSmearedFourShell(smear_coef_FourShell,smeared_FourShell,smeared_qvals)	//this wrapper fills the STRUCT
	
	Display smeared_FourShell 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("SmearedFourShell","smear_coef_FourShell","smear_parameters_FourShell","FourShell")
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 OneShell(cw,yw,xw) : FitFunc
	Wave cw,yw,xw
	
#if exists("OneShellX")
	yw = OneShellX(cw,xw)
#else
	yw = fOneShell(cw,xw)
#endif
	return(0)
End

Function TwoShell(cw,yw,xw) : FitFunc
	Wave cw,yw,xw
	
#if exists("TwoShellX")
	yw = TwoShellX(cw,xw)
#else
	yw = fTwoShell(cw,xw)
#endif
	return(0)
End

Function ThreeShell(cw,yw,xw) : FitFunc
	Wave cw,yw,xw
	
#if exists("ThreeShellX")
	yw =ThreeShellX(cw,xw)
#else
	yw = fThreeShell(cw,xw)
#endif
	return(0)
End

Function FourShell(cw,yw,xw) : FitFunc
	Wave cw,yw,xw
	
#if exists("FourShellX")
	yw = FourShellX(cw,xw)
#else
	yw = fFourShell(cw,xw)
#endif
	return(0)
End


//
// unsmeared model calculation
//
Function fOneShell(w,x) : FitFunc
	Wave w
	Variable x
	
	// variables are:
	//[0] scale factor
	//[1] radius of core [‐
	//[2] SLD of the core	[�-2]
	//[3] thickness of the shell	[‐
	//[4] SLD of the shell
	//[5] SLD of the solvent
	//[6] background	[cm-1]
	
	// All inputs are in ANGSTROMS
	//OUTPUT is normalized by the particle volume, and converted to [cm-1]
	
	
	Variable scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg
	scale = w[0]
	rcore = w[1]
	rhocore = w[2]
	thick = w[3]
	rhoshel = w[4]
	rhosolv = w[5]
	bkg = w[6]
	
	// calculates scale *( f^2 + bkg)
	Variable bes,f,vol,qr,contr,f2
	
	// core first, then add in shell
	qr=x*rcore
	contr = rhocore-rhoshel
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*rcore^3
	f = vol*bes*contr
	//now the shell
	qr=x*(rcore+thick)
	contr = rhoshel-rhosolv
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*(rcore+thick)^3
	f += vol*bes*contr
	
	// normalize to particle volume and rescale from [�-1] to [cm-1]
	f2 = f*f/vol*1.0e8
	
	//scale if desired
	f2 *= scale
	// then add in the background
	f2 += bkg
	
	return (f2)
End


Function fTwoShell(w,x) : FitFunc
	Wave w
	Variable x
	
	// variables are:
	//[0] scale factor
	//[1] radius of core [‐
	//[2] SLD of the core	[�-2]
	//[3] thickness of shell 1 [‐
	//[4] SLD of shell 1
	//[5] thickness of shell 2 [‐
	//[6] SLD of shell 2
	//[7] SLD of the solvent
	//[8] background	[cm-1]
	
	// All inputs are in ANGSTROMS
	//OUTPUT is normalized by the particle volume, and converted to [cm-1]
	
	
	Variable scale,rcore,thick1,thick2,rhocore,rhoshel1,rhoshel2,rhosolv,bkg
	scale = w[0]
	rcore = w[1]
	rhocore = w[2]
	thick1 = w[3]
	rhoshel1 = w[4]
	thick2 = w[5]
	rhoshel2 = w[6]
	rhosolv = w[7]
	bkg = w[8]
	
	// calculates scale *( f^2 + bkg)
	Variable bes,f,vol,qr,contr,f2
	
	// core first, then add in shells
	qr=x*rcore
	contr = rhocore-rhoshel1
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*rcore^3
	f = vol*bes*contr
	//now the shell (1)
	qr=x*(rcore+thick1)
	contr = rhoshel1-rhoshel2
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*(rcore+thick1)^3
	f += vol*bes*contr
	//now the shell (2)
	qr=x*(rcore+thick1+thick2)
	contr = rhoshel2-rhosolv
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*(rcore+thick1+thick2)^3
	f += vol*bes*contr
	
	// normalize to particle volume and rescale from [�-1] to [cm-1]
	f2 = f*f/vol*1.0e8
	
	//scale if desired
	f2 *= scale
	// then add in the background
	f2 += bkg
	
	return (f2)
	
End

Function fThreeShell(w,x) : FitFunc
	Wave w
	Variable x
	
	// variables are:
	//[0] scale factor
	//[1] radius of core [‐
	//[2] SLD of the core	[�-2]
	//[3] thickness of shell 1 [‐
	//[4] SLD of shell 1
	//[5] thickness of shell 2 [‐
	//[6] SLD of shell 2
	//[7] SLD of the solvent
	//[8] background	[cm-1]
	
	// All inputs are in ANGSTROMS
	//OUTPUT is normalized by the particle volume, and converted to [cm-1]
	
	
	Variable scale,rcore,thick1,thick2,thick3,rhoshel1,rhoshel2,rhoshel3
	Variable rhocore,rhosolv,bkg
	scale = w[0]
	rcore = w[1]
	rhocore = w[2]
	thick1 = w[3]
	rhoshel1 = w[4]
	thick2 = w[5]
	rhoshel2 = w[6]
	thick3 = w[7]
	rhoshel3 = w[8]
	rhosolv = w[9]
	bkg = w[10]
	
	// calculates scale *( f^2 + bkg)
	Variable bes,f,vol,qr,contr,f2
	
	// core first, then add in shells
	qr=x*rcore
	contr = rhocore-rhoshel1
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*rcore^3
	f = vol*bes*contr
	//now the shell (1)
	qr=x*(rcore+thick1)
	contr = rhoshel1-rhoshel2
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*(rcore+thick1)^3
	f += vol*bes*contr
	//now the shell (2)
	qr=x*(rcore+thick1+thick2)
	contr = rhoshel2-rhoshel3
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*(rcore+thick1+thick2)^3
	f += vol*bes*contr
	//now the shell (3)
	qr=x*(rcore+thick1+thick2+thick3)
	contr = rhoshel3-rhosolv
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*(rcore+thick1+thick2+thick3)^3
	f += vol*bes*contr
	
	// normalize to particle volume and rescale from [�-1] to [cm-1]
	f2 = f*f/vol*1.0e8
	
	//scale if desired
	f2 *= scale
	// then add in the background
	f2 += bkg
	
	return (f2)
End


Function fFourShell(w,x) : FitFunc
	Wave w
	Variable x
	
	// variables are:
	//[0] scale factor
	//[1] radius of core [‐
	//[2] SLD of the core	[�-2]
	//[3] thickness of shell 1 [‐
	//[4] SLD of shell 1
	//[5] thickness of shell 2 [‐
	//[6] SLD of shell 2
	//[7] SLD of the solvent
	//[8] background	[cm-1]
	
	// All inputs are in ANGSTROMS
	//OUTPUT is normalized by the particle volume, and converted to [cm-1]
	
	
	Variable scale,rcore,thick1,thick2,thick3,thick4
	Variable rhoshel1,rhoshel2,rhoshel3,rhoshel4
	Variable rhocore,rhosolv,bkg
	scale = w[0]
	rcore = w[1]
	rhocore = w[2]
	thick1 = w[3]
	rhoshel1 = w[4]
	thick2 = w[5]
	rhoshel2 = w[6]
	thick3 = w[7]
	rhoshel3 = w[8]
	thick4 = w[9]
	rhoshel4 = w[10]
	rhosolv = w[11]
	bkg = w[12]
	
	// calculates scale *( f^2 + bkg)
	Variable bes,f,vol,qr,contr,f2
	
	// core first, then add in shells
	qr=x*rcore
	contr = rhocore-rhoshel1
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*rcore^3
	f = vol*bes*contr
	//now the shell (1)
	qr=x*(rcore+thick1)
	contr = rhoshel1-rhoshel2
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*(rcore+thick1)^3
	f += vol*bes*contr
	//now the shell (2)
	qr=x*(rcore+thick1+thick2)
	contr = rhoshel2-rhoshel3
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*(rcore+thick1+thick2)^3
	f += vol*bes*contr
	//now the shell (3)
	qr=x*(rcore+thick1+thick2+thick3)
	contr = rhoshel3-rhoshel4
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*(rcore+thick1+thick2+thick3)^3
	f += vol*bes*contr
	//now the shell (4)
	qr=x*(rcore+thick1+thick2+thick3+thick4)
	contr = rhoshel4-rhosolv
	bes = 3*(sin(qr)-qr*cos(qr))/qr^3
	vol = 4*pi/3*(rcore+thick1+thick2+thick3+thick4)^3
	f += vol*bes*contr
	
	
	// normalize to particle volume and rescale from [�-1] to [cm-1]
	f2 = f*f/vol*1.0e8
	
	//scale if desired
	f2 *= scale
	// then add in the background
	f2 += bkg
	
	return (f2)
End




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

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

	return(0)
End

Function SmearedTwoShell(s) : FitFunc
	Struct ResSmearAAOStruct &s

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

	return(0)
End

Function SmearedThreeShell(s) : FitFunc
	Struct ResSmearAAOStruct &s

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

	return(0)
End

Function SmearedFourShell(s) : FitFunc
	Struct ResSmearAAOStruct &s

//	the name of your unsmeared model (AAO) is the first argument
	Smear_Model_20(FourShell,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 fSmearedOneShell(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 = SmearedOneShell(fs)
	
	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 fSmearedTwoShell(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 = SmearedTwoShell(fs)
	
	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 fSmearedThreeShell(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 = SmearedThreeShell(fs)
	
	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 fSmearedFourShell(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 = SmearedFourShell(fs)
	
	return (0)
End