#pragma rtGlobals=1		// Use modern global access method.
#pragma version = 1.0



// Procedures for convolution and convolution based line broadening with awareness of X-scaling in waveform data
//
// The procedures are based on the built in function "Convolve" which is FFT based, and thus fast for large data sets, but is 
// ignorant of waveform x-scaling.  Meaninful results from use of Convolve with X-scaling awareness require that any two input waves 
// (srcWaveName (kernal) and DestWaveName in the "Convolve" documentation) have the same delta X property.  The x_0 value also needs to be set correctly
// depending on the nature of the broadening and which Convolve flags are used.  The other important consideration handled by these procedures is 
// normalization.  If you are familiar with "Convolution" as a calculus concept, it is usually described as a single integral over the product of two functions, 
// and if a function is broadened with a kernal with unit area (integral = 1) then the output will preserve its normalization of the destination.  
// This is how we spectroscopists typically conceptualize "Convolve".  Because the integral has an x-spacing component to it which Igor's built-in Convolve ignores,
// normalization in the "desired sense" (meaning area returned from an integration) is preserved by scaling the kernal so its SUM over all elements is 1, 
// rather than its area (integration).   
//
// This procedure file also supplies a conversion function called "FillDestWave" (and wrapper procedure) for converting X,Y pairs to waveforms.  This conversion differs from 
// the wavemetrics supplied package "XY to waveform" in the assumption of the nature of the X,Y data.  Igor's built in "XY to waveform" assumes that the X,Y
// pairs reflect a smooth function of Y vs X with points sampled at uneven intervals, and the application is just re-sampling to a constant X spacing so
// the data can be displayed as a waveform.  If you would typically plot the X,Y data as a "lines between points" and the converted waveform produces an
// identical or very close recreation then the built in XY to Waveform is an appropriate conversion.  FillDestWave assumes the input X,Y data represents a series of
// delta functions: infinitely narrow peaks with a finite area underneath.  This comes up in spectroscopy simulations frequently because calculation of line strength and
// position is often done first and simple to apply before application of line shape.  I typically plot this data type as "Sticks to Zero" as it best visually 
// represents the data.  The difference in output is readily apparent wherever there are gaps in the X data significantly larger than the output waveform spacing.  The FillDestWave
// is similar to a weighted histogram (although there is bleeding of intensity into the nearest bin to better handle X values that are not perfectly centered in a bin)
// , where line data are summed into small(ish) bins so that the new waveform has an area integral from x = a to b matching the sum of Y values from the X,Y data for points a < X < b.
// 
// The functions supplied (FillDestWave and <shape>convolve) can all be run multiple times on the same output wave with meaningful results:  
//  The FillDestWave function can be called on the same
//  destination wave and the results will be added to existing content (makes it possible to include multiple pairs of X,Y waves into a single combined destination wave
// without needing to run a concatentation).  Calling a <shape>convolve function on an already broadened function gives the results of cumulative broadenings.  I.E., if 
// the starting data is a scan from a high resolution instrument that is already stored as a waveform, and you want to simulate a lower resolution, the additional convolution will work.
//
//  The <shape>Broaden functions, for which this procedure was originally worked out combing FillDestWave and <shape>convolve to produce line broadened waveform data from X,Y inputs
//
//  Version 1.0
//  Ian Konen
//  iankonen@gmail.com
//
static constant DEFPPHW = 5	                       // Default Points Per Half Width for setting the output wave X-scaling
static constant DEFGAUSSHWS = 5		   //  Default Halfwidths for gaussian broadening
static constant DEFLORHWS = 1000              // Default halfwidths for Lorentzian broadening
static constant GSIG2HWHM = 1.177410023   //  convert gaussian "sigma" parameter (standard deviation) to a halfwidth
static constant MAXPOINTSCONV = 1e7	   // when automatically calculating point spacing over a fixed range, don't make waves longer than this number

Menu "Line Broadening"
	 	"Gauss Convolve", GaussConvolveWrapper()
		"Lorentzian Convolve", LorConvolveWrapper()
		"Voigt Convolve", VoigtConvolveWrapper()
		"Fill Destination", FillDestWrapper()
		"Gauss Broaden", GaussBroadenWrapper()
		"Lorentzian Broaden", LorBroadenWrapper()
		"Voigt Broaden", VoigtBroadenWrapper()
End

//  Wrapper procedures for menu selection:
Proc GaussConvolveWrapper(sourcewavenm, destinationwavenm, HWHM, downsample)
	string sourcewavenm
	prompt sourcewavenm, "Source Wave", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string DestinationWavenm
	prompt DestinationWavenm, "Destination (empty overwrites source)"
	variable HWHM
	prompt HWHM, "Half Width at Half Max"
	variable downsample = 1
	prompt downsample, "Down Sampling (1 for no change)"
	if (strlen(DestinationWavenm))
		GaussConvolve($sourcewavenm,hwhm,outwavename = DestinationWavenm, downsample = downsample)
	else
		GaussConvolve($sourcewavenm,hwhm, downsample = downsample)
	endif
EndMacro

Proc LorConvolveWrapper(sourcewavenm, destinationwavenm, HWHM, downsample)
	string sourcewavenm
	prompt sourcewavenm, "Source Wave", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string DestinationWavenm
	prompt DestinationWavenm, "Destination (empty overwrites source)"
	variable HWHM
	prompt HWHM, "Half Width at Half Max"
	variable downsample = 1
	prompt downsample, "Down Sampling (1 for no change)"
	if (strlen(DestinationWavenm))
		LorConvolve($sourcewavenm,hwhm,outwavename = DestinationWavenm, downsample = downsample)
	else
		LorConvolve($sourcewavenm,hwhm, downsample = downsample)
	endif
EndMacro

Proc VoigtConvolveWrapper(sourcewavenm, destinationwavenm, L_HWHM, G_HWHM, downsample)
	string sourcewavenm
	prompt sourcewavenm, "Source Wave", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string DestinationWavenm
	prompt DestinationWavenm, "Destination (empty overwrites source)"
	variable L_HWHM
	prompt L_HWHM, "Lorentzian Width (HWHM)"
	variable G_HWHM
	prompt G_HWHM, "Guassian Width (HWHM)"
	variable downsample = 1
	prompt downsample, "Down Sampling (1 for no change)"
	if (strlen(DestinationWavenm))
		VoigtConvolve($sourcewavenm,L_hwhm,G_hwhm,outwavename = DestinationWavenm, downsample = downsample)
	else
		VoigtConvolve($sourcewavenm,L_hwhm,G_hwhm, downsample = downsample)
	endif
End

Proc GaussBroadenWrapper(Xdata, Ydata, outwavenm, HWHM)
	string Xdata
	prompt Xdata, "Wave of Line X values", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string Ydata
	prompt Ydata, "Wave of Line Y values", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string outwavenm = ""
	prompt outwavenm, "Name of Output Wave"
	variable HWHM
	prompt HWHM, "Width (HWHM)"
	GaussBroaden($Xdata, $Ydata, outwavenm, HWHM)
End

Proc LorBroadenWrapper(Xdata, Ydata, outwavenm, HWHM)
	string Xdata
	prompt Xdata, "Wave of Line X values", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string Ydata
	prompt Ydata, "Wave of Line Y values", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string outwavenm = ""
	prompt outwavenm, "Name of Output Wave"
	variable HWHM
	prompt HWHM, "Width (HWHM)"
	LorBroaden($Xdata, $Ydata, outwavenm, HWHM)
End

Proc VoigtBroadenWrapper(Xdata, Ydata, outwavenm, L_HWHM, G_HWHM)
	string Xdata
	prompt Xdata, "Wave of Line X values", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string Ydata
	prompt Ydata, "Wave of Line Y values", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string outwavenm = ""
	prompt outwavenm, "Name of Output Wave"
	variable L_HWHM
	prompt L_HWHM, "Lorentzian Width (HWHM)"
	variable G_HWHM
	prompt G_HWHM, "Gaussian Width (HWHM)"
	VoigtBroaden($Xdata, $Ydata, outwavenm, L_HWHM, G_HWHM)
End

Proc FillDestWrapper(sourceX, sourceY, destinationwavenm, xspacing, xextension)
	string sourceX
	prompt sourceX, "X data", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string sourceY
	prompt sourceY, "Y data", popup, WaveList("*",";","DIMS:1,TEXT:0")
	string destinationwavenm
	prompt destinationwavenm, "Destination Wave"
	variable xspacing = 1
	prompt xspacing, "Destination X Spacing"
	variable xextension = 1
	prompt xextension, "Extend Destination X scale by "
	variable destminx = wavemin($sourceX) - xextension
	variable destmaxx = wavemax($sourceX) + xextension
	variable destnpnts = (destmaxx - destminx) / xspacing + 1
	make /d /o /n = (destnpnts) $destinationwavenm
	SetScale /P x, (destminx), (xspacing), $destinationwavenm
	$destinationwavenm = 0
	FillDestWave($sourceX, $sourceY, $destinationwavenm)
End
 
// END Wrapper Functions
//


//  CONVOLVE FUNCTIONS

// In the following functions, the built-in Igor operation "convolve" is used
// to, well, convolve "startwave" with a peak function kernel (I've coded Gaussian, Lorentzian and Voight, but the principal can 
// be applied to any peak function.  These functions just coerce the scaling to be appropriate for waveform data (the 
// built-in Convolve operation ignores x-scaling and does everything point-to-point) and fix the normalization (to preserve Area normalization).
//
// Startwave could be an initial waveform from FillDestWave, or an already broadened wave (input an "intrinsic" signal and
// apply an instrument resolution function, for example).
//
// If the optional string "outwavename" is not supplied, the startwave is overwritten.  
// Extendhws is a measure of how may half-widths the kernel should be.  The smaller the kernel, the faster and less memory 
// intensive is the broadening operation, but tail-effects will truncated.  The default values I supplied in the constants tuncate
// the tails at < 1% maximum.  
//
//  Use downsample to decimate the wave (default behavior preserves points)
//   A negative downsample is applied BEFORE convolution (might be slightly less accurate, but better memory conservation)
//   Should use built-in low pass filtering of Resample to preserve normalization and not miss peaks.


// Convolves an existing waveform with a [further] Gaussian function.  
Function GaussConvolve(startwave,hwhm,[extendhws,outwavename,downsample])
	wave startwave
	variable hwhm
	variable extendhws
	string outwavename
	variable downsample
	if (paramisdefault(extendhws))
		extendhws = DEFGAUSSHWS
	endif
	
	if (paramisdefault(outwavename))
		wave outwave = startwave
	else
		duplicate /o startwave, $outwavename
		wave outwave = $outwavename
	endif	
	
	if (!paramisdefault(downsample))
		if (downsample < 0)			// downsample before convolution
			Resample /DOWN=(-downsample)  outwave
		endif
	endif

	variable hwpnts = ceil(extendhws*hwhm / deltax(outwave))	// halfwidth in integer data points
	make /FREE /d /n = (2*hwpnts + 1) tempgauss
	setscale /p x, (-hwpnts*deltax(outwave)), (deltax(outwave)), "", tempgauss

	tempgauss = GSIG2HWHM/sqrt(pi)/ hwhm*deltax(outwave) * exp(-(x * GSIG2HWHM / hwhm)^2)
	
	Convolve /A tempgauss, outwave
	if (!paramisdefault(downsample))
		if (downsample > 0)			// downsample after convolution
 			Resample/DOWN=(downsample) outwave
		endif
	endif

	return 1
End

// Convolves an existing waveform with a [further] Lorentzian function.  
Function LorConvolve(startwave,hwhm,[extendhws,outwavename,downsample])
	wave startwave
	variable hwhm
	variable extendhws
	string outwavename
	variable downsample
	if (paramisdefault(extendhws))
		extendhws = DEFLORHWS
	endif
	if (paramisdefault(outwavename))
		wave outwave = startwave
	else
		duplicate /o startwave, $outwavename
		wave outwave = $outwavename
		SetFormula outwave, ""
	endif
	
	if (!paramisdefault(downsample))
		if (downsample < 0)			// downsample before convolution
			Resample /DOWN=(-downsample)  outwave
		endif
	endif

	variable hwpnts = ceil(extendhws*hwhm / deltax(outwave))	// halfwidth in integer data points
	if (hwpnts > MAXPOINTSCONV)
		DoAlert 0, "Convolution requires " + num2str(2*hwpnts + 1) + "points."
		return 0
	endif
	make /FREE /d /n = (2*hwpnts + 1) templor
	setscale /p x, (-hwpnts*deltax(outwave)), (deltax(outwave)), "", templor
	templor = 1/(x^2 + hwhm^2)
	variable thesum = sum(templor)
	templor /= thesum
	Convolve /A templor, outwave
	if (!paramisdefault(downsample))
		if (downsample > 0)			// downsample before convolution
			
 			Resample/DOWN=(downsample) outwave
		endif
	endif
	return 1
End

// Actually a combined convolution of Lorentzian and Guassian (which is where the Voight function comes from).  
// Note that the input shape parameters are Gaussian width and Lorentzian width as opposed to the separate width and shape 
// parameters used by VoigtFit (from the multipeak fitting XOP).  

Function VoigtConvolve(startwave,L_hwhm,G_hwhm,[Lextendhws,Gextendhws,outwavename,downsample])
	wave startwave
	variable L_hwhm,G_hwhm
	variable Lextendhws,Gextendhws
	string outwavename
	variable downsample
	
	if (paramisdefault(Lextendhws))
		Lextendhws = DEFLORHWS
	endif
	if (paramisdefault(Gextendhws))
		Gextendhws = DEFGAUSSHWS
	endif

	if (paramisdefault(outwavename))
		wave outwave = startwave
	else
		duplicate /o startwave, $outwavename
		wave outwave = $outwavename
	endif
	
	if (!paramisdefault(downsample))
		if (downsample < 0)			// downsample before convolution
			Resample /DOWN=(-downsample) outwave
		endif
	endif

	variable hwpnts = ceil(Lextendhws*L_hwhm / deltax(startwave))	// halfwidth in integer data points, rounded up.  
	make /FREE /d /n = (2*hwpnts + 1) templor
	setscale /p x, (-hwpnts*deltax(startwave)), (deltax(startwave)), "", templor	
	
	templor = 1/(x^2 + L_hwhm^2)

	hwpnts = ceil(Gextendhws*G_hwhm / deltax(startwave))	// halfwidth in integer data points
	make /FREE /d /n = (2*hwpnts + 1) tempgauss
	setscale /p x, (-hwpnts*deltax(startwave)), (deltax(startwave)), "", tempgauss	

	tempgauss = exp(-(x * GSIG2HWHM / G_hwhm)^2)
	Convolve tempgauss,templor		// templor is now really a tempVoight
	variable thesum = sum(templor)	
	templor /= thesum				// Now it's a normalization preserving Voight
	Convolve /A templor, outwave
	
	if (!paramisdefault(downsample))
		if (downsample > 0)			// downsample after convolution			
 			Resample/DOWN=(downsample) outwave
		endif
	endif

	return 1
End

// END CONVOLVE FUNCTIONS



// BROADENING FUNCTIONS


// NOTE : destination wave is not pre-cleared by setting it to 0.  This should typically be done by the programmer or user
// previous to calling this function, but it is possible to call this multiple times using
// different X,Y pairs to build up the contents of a dest wave that represents a total signature.

Function FillDestWave(pos, inten, dest)
	wave pos, inten, dest
	variable cnt = 0, plow, lowint, highint
	variable leftex = leftx(dest)
	variable deltaex = deltax(dest)
	variable endex = rightx(dest) - deltaex
//	dest = 0
	variable loopmax = min(numpnts(pos),numpnts(inten))
	for (cnt = 0; cnt < loopmax; cnt += 1)
		if (!(pos[cnt] < leftex || pos[cnt] > endex))
			plow = (pos[cnt] - leftx(dest))/deltaex
//			highint = inten[cnt]*(plow - floor(plow))
//			lowint = inten[cnt] - highint
			highint = inten[cnt]*(plow - floor(plow))/deltaex
			lowint = inten[cnt]/deltaex - highint
			plow = floor(plow)
			dest[plow] += lowint
			dest[plow + 1] += highint
		endif
	endfor
End


Function GaussBroaden(Xwave, Ywave, outwavenm, hwhm, [LowX, HighX, Xspacing])
	wave Xwave, Ywave
	string outwavenm
	variable hwhm
	variable LowX, HighX, Xspacing
	if (ParamIsDefault(LowX))
		LowX = WaveMin(Xwave) - hwhm*5
	endif
	if (ParamIsDefault(HighX))
		HighX = WaveMax(Xwave) + hwhm*5
	endif
	if (ParamIsDefault(Xspacing))
		Xspacing = hwhm / DEFPPHW
	endif
	variable numpoints = (HighX - LowX) / Xspacing + 1
	if (numpoints > MAXPOINTSCONV)
		print "Output Wave requires ", numpoints, " points."
		print "Override defaults to make a smaller wave, or manually fill and broaden."
		return 0
	endif
	make /o /d /n = (numpoints) $outwavenm
	wave outwave = $outwavenm
	outwave = 0
	SetScale /P x, (LowX), (Xspacing), "", outwave
	FillDestWave(Xwave, Ywave, outwave)
	GaussConvolve(outwave,hwhm)
	return 1
End

Function LorBroaden(Xwave, Ywave, outwavenm, hwhm, [LowX, HighX, Xspacing])
	wave Xwave, Ywave
	string outwavenm
	variable hwhm
	variable LowX, HighX, Xspacing
	if (ParamIsDefault(LowX))
		LowX = WaveMin(Xwave) - hwhm*5
	endif
	if (ParamIsDefault(HighX))
		HighX = WaveMax(Xwave) + hwhm*5
	endif
	if (ParamIsDefault(Xspacing))
		Xspacing = hwhm / DEFPPHW
	endif
	variable numpoints = (HighX - LowX) / Xspacing + 1
	if (numpoints > MAXPOINTSCONV)
		print "Default Output results in ", numpoints, " points."
		print "Override defaults to make a smaller wave, or manually fill and broaden."
		return 0
	endif
	make /o /d /n = (numpoints) $outwavenm
	wave outwave = $outwavenm
	outwave = 0
	SetScale /P x, (LowX), (Xspacing), "", outwave
	FillDestWave(Xwave, Ywave, outwave)
	LorConvolve(outwave,hwhm)
	return 1
End

Function VoigtBroaden(Xwave, Ywave, outwavenm, l_hwhm, g_hwhm, [LowX, HighX, Xspacing])
	wave Xwave, Ywave
	string outwavenm
	variable l_hwhm, g_hwhm
	variable LowX, HighX, Xspacing
	variable v_hwhm = 0.5*l_hwhm + sqrt(0.25*l_hwhm^2 + g_hwhm^2)	// hardly need this, could just use a minmax argument
	if (ParamIsDefault(LowX))
		LowX = WaveMin(Xwave) - v_hwhm*5
	endif
	if (ParamIsDefault(HighX))
		HighX = WaveMax(Xwave) + v_hwhm*5
	endif
	if (ParamIsDefault(Xspacing))
		Xspacing = v_hwhm / DEFPPHW
	endif
	variable numpoints = (HighX - LowX) / Xspacing + 1
	if (numpoints > MAXPOINTSCONV)
		print "Default Output results in ", numpoints, " points."
		print "Override defaults to make a smaller wave, or manually fill and broaden."
		return 0
	endif
	make /o /d /n = (numpoints) $outwavenm
	wave outwave = $outwavenm
	outwave = 0
	SetScale /P x, (LowX), (Xspacing), "", outwave
	FillDestWave(Xwave, Ywave, outwave)
	VoigtConvolve(outwave,l_hwhm, g_hwhm)
	return 1
End