#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3
#pragma version=1.00

// for peak fitting, a collection of lineshapes 
// and wrapper functions for fitting

// https://www.wavemetrics.com/user/tony

static constant ln2=0.69314718055994528623
static constant sqrtln2=0.832554611157698			// sqrt(ln(2))
static constant sqrtln2pi=0.469718639349826		// sqrt(ln(2)/pi)


// Gaussian-Lorentzian Cross Product (Amplitude)
// a0 = amplitude
// a1
 = center
// a2 = width (>0) 
// a3 = shape (≥0, ≤1)
function GL_CrossProduct(x, a0, a1, a2, a3)
	variable x, a0, a1, a2, a3
		
	return a0/(1+a3*((x-a1)/a2)^2*exp((1-a3)/2*((x-a1)/a2)^2))	
end


// Gaussian-Lorentzian Sum (Amplitude)
// a0= amplitude
// a1
= center
// a2= width (>0)
// a3 = shape (≥0, ≤1)
function GL_Sum_amp(x, a0, a1, a2, a3)
	variable x, a0, a1, a2, a3
		
	return a0*( a3/a2*sqrtln2pi*exp(-4*ln2*((x-a1)/a2)^2) + ((1-a3)/(pi*a2*(1+4*((x-a1)/a2)^2))) / (a3/a2*sqrtln2pi + (1-a3)/(pi*a2)) )
end


// Gaussian-Lorentzian Sum (Area)
// a0= area
// a1
= center
// a2= width (>0) 
// a3 = shape (≥0, ≤1)
function GL_Sum_area(x, a0, a1, a2, a3)
	variable x, a0, a1, a2, a3
		
	return 2*a0*( a3/a2*sqrtln2pi*exp(-4*ln2*((x-a1)/a2)^2)  +  (1-a3)/(pi*a2*(1+4*((x-a1)/a2)^2)) )
end


// A Voigt approximation
// Wertheim, G.K., Butler, M.A., West, K.W., and Buchanan, D.N.E. (1974).
// Determination of the Gaussian and Lorentzian content of experimental 
// line shapes. Review of Scientific Instruments. 45 (11): 1369–1371. 
function Wertheim(G, nu, x, x0)
	variable G, nu, x, x0
	// G is width of pseudo-Voigt
	
	variable b=sqrtln2/2
	
	return nu/(1+((x-x0)/(0.5*G))^2) + (1-nu)*exp(-((x-x0)/(b*G))^2)	
end

// Two Voigt approximations from CasaXPS
// m is mixing term
// E is centre
// http://www.casaxps.com/help_manual/line_shapes.htm
// I think these are the same as the cross-product and sum functions with some of the constants wrapped up 
function CASA_GL(x, F, E, m)
	variable x, F, E, m
	
	return exp(-4*ln(2)*(1-m)*(x-E)^2/F^2)/(1+4*m*(x-E)^2/F^2)
end

function CASA_SGL(x, F, E, m)
	variable x, F, E, m
	
	return (1-m)*exp(-4*ln(2)*(x-E)^2/F^2) + m/(1+4*(x-E)^2/F^2)
end


// -------------- wrappers for these functions for curve fit dialog --------------

Function VoigtFit(w,xx) : FitFunc
	Wave w
	Variable xx

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(xx) = Y0+Amp*VoigtFunc(width*(xx-x0),shape)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ xx
	//CurveFitDialog/ Coefficients 5
	//CurveFitDialog/ w[0] = Y0
	//CurveFitDialog/ w[1] = Amp
	//CurveFitDialog/ w[2] = width
	//CurveFitDialog/ w[3] = x0
	//CurveFitDialog/ w[4] = shape
	
//	a = w_coef[1]*sqrt(pi)/w_coef[2]
//	wg = sqrt(ln(2))/w_coef[2]
//	wl = w_coef[4]/w_coef[2] 
//	wv = wl/2 + sqrt( wl^2/4 + wg^2)

//	Wells, R.J., Rapid Approximation to the Voigt/Faddeeva Function 
// and Its Derivatives, Journal of Quantitative Sprectroscopy and 
// Radiative Transfer, 1999

	return w[0]+w[1]*VoigtFunc(w[2]*(xx-w[3]),w[4])
End

Function VoigtFitNoBase(w,xx) : FitFunc
	Wave w
	Variable xx

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(xx) = Amp*VoigtFunc(width*(xx-x0),shape)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ xx
	//CurveFitDialog/ Coefficients 4
	//CurveFitDialog/ w[0] = Amp
	//CurveFitDialog/ w[1] = width
	//CurveFitDialog/ w[2] = x0
	//CurveFitDialog/ w[3] = shape
	
//	a = w_coef[0]*sqrt(pi)/w_coef[1]
//	wg = sqrt(ln(2))/w_coef[1]
//	wl = w_coef[3]/w_coef[1] 
//	wv = wl/2 + sqrt( wl^2/4 + wg^2)
	
	return w[0]*VoigtFunc(w[1]*(xx-w[2]),w[3])
End



Function fitCasaSGL(w,x) : FitFunc
	Wave w
	Variable x

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(x) = amp*CASA_SGL(x,width,centre,m)+a+b*x
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 6
	//CurveFitDialog/ w[0] = width
	//CurveFitDialog/ w[1] = centre
	//CurveFitDialog/ w[2] = m
	//CurveFitDialog/ w[3] = amp
	//CurveFitDialog/ w[4] = a
	//CurveFitDialog/ w[5] = b

	return w[3]*CASA_SGL(x,w[0],w[1],w[2])+w[4]+w[5]*x
End

Function fitCasaGL(w,x) : FitFunc
	Wave w
	Variable x

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(x) = amp*CASA_GL(x,width,centre,m)+a+b*x
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 6
	//CurveFitDialog/ w[0] = width
	//CurveFitDialog/ w[1] = centre
	//CurveFitDialog/ w[2] = m
	//CurveFitDialog/ w[3] = amp
	//CurveFitDialog/ w[4] = a
	//CurveFitDialog/ w[5] = b

	return w[3]*CASA_GL(x,w[0],w[1],w[2])+w[4]+w[5]*x
End


Function fitWertheim(w,x) : FitFunc
	Wave w
	Variable x

	//CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
	//CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
	//CurveFitDialog/ Equation:
	//CurveFitDialog/ f(x) = amp*Wertheim(G, nu, x, x0) + a +b*x
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 6
	//CurveFitDialog/ w[0] = amp
	//CurveFitDialog/ w[1] = G
	//CurveFitDialog/ w[2] = nu
	//CurveFitDialog/ w[3] = x0
	//CurveFitDialog/ w[4] = a
	//CurveFitDialog/ w[5] = b

	return w[0]*Wertheim(w[1], w[2], x, w[3]) + w[4] +w[5]*x
End