#pragma rtGlobals=1		// Use modern global access method.
static constant PeakFuncInfo_ParamNames			= 0 														 // keep a separate instance of the switch constants
static constant PeakFuncInfo_PeakFName			= 1														 // allows to load this procedure separately
static constant PeakFuncInfo_GaussConvFName		= 2
static constant PeakFuncInfo_ParameterFunc		= 3
static constant PeakFuncInfo_DerivedParamNames	= 4

Function/S DS_PeakFuncInfo(InfoDesired)
	Variable InfoDesired
 
	String info=""
 
	Switch (InfoDesired)
		case PeakFuncInfo_ParamNames:
			info = "Location;Width;Height;Asymmetry;GaussW;"
			break;
		case PeakFuncInfo_PeakFName:
			info = "DSPeak"
			break;
		case PeakFuncInfo_GaussConvFName:
			info = "GausstoDSGuess"
			break;
		case PeakFuncInfo_ParameterFunc:
			info = "DSParams"
			break;
		case PeakFuncInfo_DerivedParamNames:
			info = "Location;Height;Area;FWHM;Asymmetry;GaussWidth;"
			break;
		default:
			break;
	endSwitch
 
	return info
End
//################################################################################################
Function DSPeak(w,yw,xw)
	Wave w
	Wave yw, xw
 
 	Variable Resolution	= 10																			 // the sample resolution (higher = slower)
	Variable Xdelta		= abs(xw[0] - xw[NumPnts(xw)-1]) / (Resolution * NumPnts(xw))						 // calculate an high-res xdelta for scaling
	Variable yPoints		= Resolution * NumPnts(yw)														 // make the wave X times bigger than the original y wave
	Make/FREE/D/N=(2*yPoints+1) yWave																	 // create a wave two times as big to have space before and after the peak
	SetScale/P x xw[0]-yPoints/2*Xdelta, Xdelta, yWave														 // reserve 1/2 size of space before the x wave starts (1/2 after)
 
	yWave = w[2]*cos(Pi*w[3]/2+(1-w[3])*atan((x-w[0])/w[1]))/((x-w[0])^2+w[1]^2)^((1-w[3])/2)
	//+++++++++++++++++++++++++++++++ convolve with gaussian +++++++++++++++++++++++++++++++++	
	if (w[4] != 0)																							 // convolve with gaussian if width is set
		Variable width	= w[4]/Xdelta																		 // width normalized by x scaling
		yPoints			= 10*width																		 // create enough points to fit in the full gaussian
		Make/FREE/D/N=(yPoints + 1) GaussWave															 // create a gaussian for the convolution
		GaussWave = Gauss(x,yPoints/2,width)
		Convolve/A GaussWave yWave																		 // convolve with the prepared gaussian
	endif
	yw = yWave(xw[p])																					 // write the calculated values into the resultwave
End
//################################################################################################
Function GaussToDSGuess(w)
	Wave w
	
	Variable x0		= w[0]																				 // extract the initial guesses
	Variable width	= w[1]
	Variable height	= w[2]
	Variable lwidth	= w[3]
	Variable rwidth	= w[4]
	
	Redimension/N=5 w																					 // redimension to apropriate values
	// output values are: "Location;Width;Height;Asym;GaussW;"
	w[1] = 0.8*width																						 // 80% into lorentzian width
	w[2] = 0.8*height																						 // make the peak a bit smaller
	w[3] = 5*(rwidth-lwidth)/width																			 // asymmetry depending on width difference
	w[4] = 0.3*width																						 // 30% into gaussian width
	return 0
end
//################################################################################################
Function DSParams(cw, sw, outWave)
	Wave cw, sw, outWave

	Variable Amplitude, Location, FWHM, PeakArea
	DS_FindStats(cw, "DSPeak", Amplitude, Location, PeakArea, FWHM)										 // approximately calculates all the results from the parameters

	Redimension/N=(6,2) outWave
	// output values are: "Location;Height;Area;FWHM;Asym;GaussWidth"	
	outWave[0][0] = Location																				 // position
	outWave[0][1] = NaN
	outWave[1][0] = Amplitude																				 // height
	outWave[1][1] = NaN
	outWave[2][0] = PeakArea																				 // area (probably wrong)
	outWave[2][1] = NaN
	outWave[3][0] = FWHM																				 // FWHM
	outWave[3][1] = NaN
	outWave[4][0] = cw[3]																					 // asymmetry value
	outWave[4][1] = sqrt(sw[3][3])
	outWave[5][0] = cw[4]																					 // gauss broadening width
	outWave[5][1] = sqrt(sw[3][3])
End
//###################################### Value finder helper function ###################################
Function DS_FindStats(cw, PeakFunc, Amplitude, Location, PeakArea, FWHM)									 // recalculate the FULL peak and extract all result values (approximately)
	Wave cw
	String PeakFunc
	Variable &Amplitude, &Location, &PeakArea, &FWHM
	
	Funcref DSPeak PeakFunction = $PeakFunc															 // reference the passed peak function (so that this function can be used for multiple peak shapes)
	
	Variable x0		= cw[0]
	Variable width	= cw[1] + cw[4]
	Variable PosH	= cw[2] > 0																			 // Look if height is positive
	
	Variable Size = 3000																					 // create a big wave (fine increments)
	Make/FREE/N=(Size+1) yw, xw
	xw = (p/Size-0.5) * 30 * width + x0																		 // create a broad x scale (to make sure the full peak fits in)
	PeakFunction(cw, yw, xw)																				 // recalculate the peak with the coef set
	WaveStats/Q yw																						 // find the statistics

	if (PosH)
		Amplitude	= V_max
		Location		= xw(V_maxloc)
	else
		Amplitude	= V_min
		Location		= xw(V_minloc)
	endif
	
	FindLevel/Q/R=(0,Size/2) yw, Amplitude/2
	Variable Left		= V_flag == 1 ? NaN : xw(V_LevelX)
	FindLevel/Q/R=(Size/2,Size) yw, Amplitude/2
	Variable Right	= V_flag == 1 ? NaN : xw(V_LevelX)	
	
	FWHM		= abs(Left -Right)
	PeakArea	= AreaXY(xw, yw)
	return 0
End