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

/////////////////////////////////////////////////
//  9/26/2011 - need to check the following 2 functions
//  does x scaling need to be the same since this is ignored in th convolution???
//////////////////////////////////////////////////

Function ExpPowerGauss(pw, yw, xw) : FitFunc
WAVE pw, yw, xw
// pw[0] = gaussian baseline offset
// pw[1] = exponential amplitude
// pw[2] = gaussian position
// pw[3] = gaussian width
// pw[4] = expoenential rate constant
// pw[5] = power function amplitude
// pw[6] = power function exponent

// convolution of a laser pulse (gaussian) with 
// an exponential (first order kietic process) and 
// a power function (higher order kinetic process)

// Arbitrary choice, make sure we are oversmapling laser pulse width
Variable resolutionFactor = 10

// make sure we oversample the gaussian impulse function
// should be the sortest time scale in our measurement

//Variable dT = min(1/(resolutionFactor*pw[4]), pw[3]/resolutionFactor)
Variable dT = pw[3]/resolutionFactor

// Calculate suitable number points for the gaussian. Length of 
//  wave is doubled so function can start 
// in the middle; +1 to make it odd so  starts at t=0, and t=0 is 
// exactly the middle point. That is better for the convolution.

// this gives us a 10X oversampled wave that is 10 widths long on each side
Variable nGaussWavePnts = round((resolutionFactor*pw[3])/dT)*2 


//Make/D/O/N=(nGaussWavePnts+1) gausswave

// In this version of the function, we make a y output wave ourselves, so 
// that we can control the resolution and accuracy of the calculation. It 
// also will allow us to use a wave assignment later to solve the problem 
// of variable X spacing or missing points.

// need plenty of points to avoid end effects when function doesn't decay to zero
// sensitive to fitting range
// be careful here!

Make/D/O/N=(nGaussWavePnts) gausswave

// This wave scaling is set such that the gaussian will 
// start at the middle of the wave 
setscale/P x -dT*(nGaussWavePnts/2),dT,gausswave

// add at least 30 sigmas onto end of wave to deal with end effects 
// so we can delete them later
Variable endp = dimsize(xw,0)
Variable xfinal = xw[(endp-1)] + pw[3]*30

//Variable nYPnts = 10*max(resolutionFactor*numpnts(yw), nGaussWavePnts) 
Variable nYPnts = 1*max(resolutionFactor*numpnts(yw), nGaussWavePnts) 
Make/C/D/O/N=(nYPnts) exppowerwave
Make/D/O/N=(nYPnts) yEPFWave

// Set the wave scaling of the intermediate output wave to have the resolution 
// calculated above, and to start at the first X value. 
//setscale/P x xw[0],dT/pw[4],powerwave
//setscale/P x xw[0],dT/pw[4], yPFWave
Variable EPFdeta = (xfinal-xw[0])/endp
Variable NdelPts = round(pw[3]*30/EPFdeta)

setscale/I x xw[0],xfinal, exppowerwave
setscale/I x xw[0],xfinal, yEPFWave


// fill expwave with exponential decay 
gausswave = (1/(pw[3]*sqrt(2*pi)))* exp(-0.5*((x)/pw[3])^2)

// Put a Gaussian peak into the intermediate output wave. We use our own wave 
// (yWave) because the convolution requires a wave with even spacing in X, 
// whereas we may get X values input that are not evenly spaced. 

// should these waves be defined for x - pw[2] > 0???
//	exppowerwave = pw[0]+(x>0)*pw[5]*((x-pw[2])^(-pw[6]))/(dT^(-pw[6])) 
//	yEPFWave =real(exppowerwave)+(x>=0)*pw[1]*(dT/pw[4])*exp(-pw[4]*(x-pw[2]))
	
	exppowerwave = pw[0]+((x-pw[2])>0)*pw[5]*((x-pw[2])^(-pw[6]))/(dT^(-pw[6])) 
	yEPFWave =real(exppowerwave)+((x-pw[2])>=0)*pw[1]*(dT/pw[4])*exp(-pw[4]*(x-pw[2]))

// and convolve with the exponential; NOTE /A 
Duplicate/O yEPFWave, yEFPo
convolve/A gausswave, yEPFWave

// Move appropriate values corresponding to input X data into the output Y 
// wave. We use a wave assignment involving the input X wave. This will 
// extract the appropriate values from the intermediate wave, interpolating 
// values where the intermediate wave doesnŐt have a value precisely at an X 
// value that is required by the input X wave. This wave assignment also
// solves the problem with auto-destination: The function can be called with 
// an X wave that sets any X spacing, so it doesnŐt matter what X values 
// are required.

// remove end effects
Variable StDel = DimSize(yEPFWave,0)-NdelPts
DeletePoints StDel, NdelPts, yEPFWave
yw = yEPFWave(xw[p])

End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////

Function PowerGauss(pw, yw, xw) : FitFunc
WAVE pw, yw, xw
// pw[0] = gaussian baseline offset
// pw[1] = gaussian amplitude
// pw[2] = gaussian position
// pw[3] = gaussian width
// pw[4] = power function exponent

// convolution of a laser pulse (gaussian) with 
// a power function (higher order kinetic process)

// Arbitrary choice
Variable resolutionFactor = 10

// make sure we oversample the gaussian impulse function
// should be the sortest time scale in our measurement

//Variable dT = min(1/(resolutionFactor*pw[4]), pw[3]/resolutionFactor)
Variable dT = pw[3]/resolutionFactor

// Calculate suitable number points for the gaussian. Length of 
//  wave is doubled so function can start 
// in the middle; +1 to make it odd so  starts at t=0, and t=0 is 
// exactly the middle point. That is better for the convolution.

Variable nGaussWavePnts = round((resolutionFactor*pw[3])/dT)*2 

//Make/D/O/N=(nGaussWavePnts+1) gausswave

// In this version of the function, we make a y output wave ourselves, so 
// that we can control the resolution and accuracy of the calculation. It 
// also will allow us to use a wave assignment later to solve the problem 
// of variable X spacing or missing points.

// need plenty of points to avoid end effects when function doesn't decay to zero
// sensitive to fitting range
// be careful here!

Make/D/O/N=(nGaussWavePnts) gausswave

// add at least 30 sigmas onto end of wave to deal with end effects 
// so we can delete them later
Variable endp = dimsize(xw,0)
Variable xfinal = xw[(endp-1)] + pw[3]*30

Variable nYPnts = 10*max(resolutionFactor*numpnts(yw), nGaussWavePnts) 
Make/C/D/O/N=(nYPnts) powerwave
Make/D/O/N=(nYPnts) yPFWave

// This wave scaling is set such that the gaussian will 
// start at the middle of the wave 
setscale/P x -dT*(nGaussWavePnts/2),dT,gausswave

// Set the wave scaling of the intermediate output wave to have the resolution 
// calculated above, and to start at the first X value. 
//setscale/P x xw[0],dT/pw[4],powerwave
//setscale/P x xw[0],dT/pw[4], yPFWave
Variable PFdeta = (xfinal-xw[0])/endp
Variable NdelPts = round(pw[3]*30/PFdeta)

setscale/I x xw[0],xfinal,powerwave
setscale/I x xw[0],xfinal, yPFWave


// fill expwave with exponential decay 
gausswave = (1/(pw[3]*sqrt(2*pi)))* exp(-0.5*((x)/pw[3])^2)

// Normalize exponential so that convolution doesn't change 
// the amplitude of the result 
//Variable sumexp 
//sumexp = sum(expwave, -inf,inf)
//expwave /= sumexp

// Put a Gaussian peak into the intermediate output wave. We use our own wave 
// (yWave) because the convolution requires a wave with even spacing in X, 
// whereas we may get X values input that are not evenly spaced. 

// same here (x-pw[2])>0
	powerwave = pw[0]+((x-pw[2])>0)*pw[1]*(x-pw[2])^(-pw[4])/(dT^(-pw[4])) 
	yPFWave = real(powerwave)

// and convolve with the exponential; NOTE /A 
convolve/A gausswave, yPFWave

// Move appropriate values corresponding to input X data into the output Y 
// wave. We use a wave assignment involving the input X wave. This will 
// extract the appropriate values from the intermediate wave, interpolating 
// values where the intermediate wave doesnŐt have a value precisely at an X 
// value that is required by the input X wave. This wave assignment also
// solves the problem with auto-destination: The function can be called with 
// an X wave that sets any X spacing, so it doesnŐt matter what X values 
// are required.

// remove end effects
Variable StDel = DimSize(yPFWave,0)-NdelPts
DeletePoints StDel, NdelPts, yPFWave
yw = yPFWave(xw[p])

End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


// Misc Fit Functions
// Keep t0 and pulsewidth identical for all decays
// G are rates

Function TwoExpGaussInf(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) = y0+(amp1/G1)*ExpGauss(x-t0,G1,pw)+(amp2/G2)*ExpGauss(x-t0,G2,pw)+(Ainf/1E-9)*ExpGauss(x-t0,1E-9,pw)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 8
	//CurveFitDialog/ w[0] = y0
	//CurveFitDialog/ w[1] = amp1
	//CurveFitDialog/ w[2] = t0
	//CurveFitDialog/ w[3] = pw
	//CurveFitDialog/ w[4] = G1
	//CurveFitDialog/ w[5] = amp2
	//CurveFitDialog/ w[6] = G2
	//CurveFitDialog/ w[7] = Ainf

	return w[0]+(w[1]/w[4])*ExpGauss(x-w[2],w[4],w[3])+(w[5]/w[6])*ExpGauss(x-w[2],w[6],w[3])+(w[7]/1E-9)*ExpGauss(x-w[2],1E-9,w[3])
End


//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


// Misc Fit Functions
// Keep t0 and pulsewidth identical for all decays
// G are rates

Function TwoExpGauss(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) = y0+(amp1/G1)*ExpGauss(x-t0,G1,pw)+(amp2/G2)*ExpGauss(x-t0,G2,pw)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 7
	//CurveFitDialog/ w[0] = y0
	//CurveFitDialog/ w[1] = amp1
	//CurveFitDialog/ w[2] = t0
	//CurveFitDialog/ w[3] = pw
	//CurveFitDialog/ w[4] = G1
	//CurveFitDialog/ w[5] = amp2
	//CurveFitDialog/ w[6] = G2

	return w[0]+(w[1]/w[4])*ExpGauss(x-w[2],w[4],w[3])+(w[5]/w[6])*ExpGauss(x-w[2],w[6],w[3])
End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function ThreeExpGaussMatt(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) = y0+(amp1/G1)*ExpGauss(x-t0,G1,pw)+(amp2/G2)*ExpGauss(x-t0,G2,pw)+(amp3/G3)*ExpGauss(x-t0,G3,pw)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 9
	//CurveFitDialog/ w[0] = y0
	//CurveFitDialog/ w[1] = amp1
	//CurveFitDialog/ w[2] = t0
	//CurveFitDialog/ w[3] = pw
	//CurveFitDialog/ w[4] = G1
	//CurveFitDialog/ w[5] = amp2
	//CurveFitDialog/ w[6] = G2
	//CurveFitDialog/ w[7] = amp3
	//CurveFitDialog/ w[8] = G3

	return w[0]+(w[1]/w[4])*ExpGauss(x-w[2],w[4],w[3])+(w[5]/w[6])*ExpGauss(x-w[2],w[6],w[3])+(w[7]/w[8])*ExpGauss(x-w[2],w[8],w[3])
End
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function FourExpGaussMatt(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) = y0+(amp1/G1)*ExpGauss(x-t0,G1,pw)+(amp2/G2)*ExpGauss(x-t0,G2,pw)+(amp3/G3)*ExpGauss(x-t0,G3,pw)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 11
	//CurveFitDialog/ w[0] = y0
	//CurveFitDialog/ w[1] = amp1
	//CurveFitDialog/ w[2] = t0
	//CurveFitDialog/ w[3] = pw	
	//CurveFitDialog/ w[4] = G1
	//CurveFitDialog/ w[5] = amp2
	//CurveFitDialog/ w[6] = G2
	//CurveFitDialog/ w[7] = amp3
	//CurveFitDialog/ w[8] = G3
	//CurveFitDialog/ w[9] = amp4
	//CurveFitDialog/ w[10] = G4
	

	return w[0]+(w[1]/w[4])*ExpGauss(x-w[2],w[4],w[3])+(w[5]/w[6])*ExpGauss(x-w[2],w[6],w[3])+(w[7]/w[8])*ExpGauss(x-w[2],w[8],w[3])+(w[9]/w[10])*ExpGauss(x-w[2],w[10],w[3])
End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function ExpGaussDampedOsc(w,t) : FitFunc
	Wave w
	Variable t

	//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(t) = (A/G1)*ExpGauss(t-t0,G1,sig)+ (B/G2)*ExpGauss(t-t0,G2,sig)+(C/G3)*ExpGauss(t-t0,G3,sig)*(D+sin(w*t+phi))
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ t
	//CurveFitDialog/ Coefficients 11
	//CurveFitDialog/ w[0] = A
	//CurveFitDialog/ w[1] = t0
	//CurveFitDialog/ w[2] = sig
	//CurveFitDialog/ w[3] = G1
	//CurveFitDialog/ w[4] = B
	//CurveFitDialog/ w[5] = G2
	//CurveFitDialog/ w[6] = C
	//CurveFitDialog/ w[7] = G3
	//CurveFitDialog/ w[8] = D
	//CurveFitDialog/ w[9] = w
	//CurveFitDialog/ w[10] = phi

	return (w[0]/w[3])*ExpGauss(t-w[1],w[3],w[2])+ (w[4]/w[5])*ExpGauss(t-w[1],w[5],w[2])+(w[6]/w[7])*ExpGauss(t-w[1],w[7],w[2])*(w[8]+sin(w[9]*t+w[10]))
End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function ExpGaussDampedOsc2(w,t) : FitFunc
	Wave w
	Variable t

	//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(t) = (A/G1)*ExpGauss(t-t0,G1,sig)+ (B/G2)*ExpGauss(t-t0,G2,sig)+(C/G3)*ExpGauss(t-t0,G3,sig)*(D+sin(w*t+phi))+(E/G4)*ExpGauss(t-t0,G4,sig)*(F+sin(w2*t+phi2))
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ t
	//CurveFitDialog/ Coefficients 16
	//CurveFitDialog/ w[0] = A
	//CurveFitDialog/ w[1] = t0
	//CurveFitDialog/ w[2] = sig
	//CurveFitDialog/ w[3] = G1
	//CurveFitDialog/ w[4] = B
	//CurveFitDialog/ w[5] = G2
	//CurveFitDialog/ w[6] = C
	//CurveFitDialog/ w[7] = G3
	//CurveFitDialog/ w[8] = D
	//CurveFitDialog/ w[9] = w
	//CurveFitDialog/ w[10] = phi
	//CurveFitDialog/ w[11] = E
	//CurveFitDialog/ w[12] = G4
	//CurveFitDialog/ w[13] = F
	//CurveFitDialog/ w[14] = w2
	//CurveFitDialog/ w[15] = phi2

	return (w[0]/w[3])*ExpGauss(t-w[1],w[3],w[2])+ (w[4]/w[5])*ExpGauss(t-w[1],w[5],w[2])+(w[6]/w[7])*ExpGauss(t-w[1],w[7],w[2])*(w[8]+sin(w[9]*t+w[10]))+(w[11]/w[12])*ExpGauss(t-w[1],w[12],w[2])*(w[13]+sin(w[14]*t+w[15]))
End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function ExpGaussDampedOsc3(w,t) : FitFunc
	Wave w
	Variable t

	//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(t) = (A/G1)*ExpGauss(t-t0,G1,sig)+ (C/G3)*ExpGauss(t-t0,G3,sig)*(D+sin(w*t+phi))
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ t
	//CurveFitDialog/ Coefficients 9
	//CurveFitDialog/ w[0] = A
	//CurveFitDialog/ w[1] = t0
	//CurveFitDialog/ w[2] = sig
	//CurveFitDialog/ w[3] = G1
	//CurveFitDialog/ w[4] = C
	//CurveFitDialog/ w[5] = G3
	//CurveFitDialog/ w[6] = D
	//CurveFitDialog/ w[7] = w
	//CurveFitDialog/ w[8] = phi

	return (w[0]/w[3])*ExpGauss(t-w[1],w[3],w[2])+ (w[4]/w[5])*ExpGauss(t-w[1],w[5],w[2])*(w[6]+sin(w[7]*t+w[8]))
End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function ExpPowerTail(w,t) : FitFunc
	Wave w
	Variable t

	//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(t) = y0+ A*exp(-t/tau)+B*t^(-alpha)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ t
	//CurveFitDialog/ Coefficients 5
	//CurveFitDialog/ w[0] = A
	//CurveFitDialog/ w[1] = tau
	//CurveFitDialog/ w[2] = B
	//CurveFitDialog/ w[3] = alpha
	//CurveFitDialog/ w[4] = y0

	return w[4]+ w[0]*exp(-t/w[1])+w[2]*t^(-w[3])
End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function Exp2PowerTail(w,t) : FitFunc
	Wave w
	Variable t

	//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(t) = y0+ A1*exp(-t/tau1)+A2*exp(-t/tau2)+B*t^(-alpha)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ t
	//CurveFitDialog/ Coefficients 7
	//CurveFitDialog/ w[0] = A1
	//CurveFitDialog/ w[1] = tau1
	//CurveFitDialog/ w[2] = A2
	//CurveFitDialog/ w[3] = tau2
	//CurveFitDialog/ w[4] = B
	//CurveFitDialog/ w[5] = alpha
	//CurveFitDialog/ w[6] = y0

	return w[6]+ w[0]*exp(-t/w[1])+w[2]*exp(-t/w[3])+w[4]*t^(-w[5])
End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function SolventResponse(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) = A* Gauss(x,t0,sig)+B*hermiteGauss(1,((x-t0)/(2*sig^2)))*Gauss(x,t0,sig)+C*hermiteGauss(2,((x-t0)/(2*sig^2)))*Gauss(x,t0,sig)
	//CurveFitDialog/ End of Equation
	//CurveFitDialog/ Independent Variables 1
	//CurveFitDialog/ x
	//CurveFitDialog/ Coefficients 5
	//CurveFitDialog/ w[0] = A
	//CurveFitDialog/ w[1] = B
	//CurveFitDialog/ w[2] = C
	//CurveFitDialog/ w[3] = t0
	//CurveFitDialog/ w[4] = sig

	return w[0]* Gauss(x,w[3],w[4])+w[1]*hermiteGauss(1,((x-w[3])/(2*w[4]^2)))*Gauss(x,w[3],w[4])+w[2]*hermiteGauss(2,((x-w[3])/(2*w[4]^2)))*Gauss(x,w[3],w[4])
End


//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function ExpConvPower(pw, yw, xw) : FitFunc 
WAVE pw, yw, xw
// pw[0] = baseline offset
// pw[1] = timezero
// pw[2] = amplitude 
// pw[3] = exponential rate constant
// pw[4] = power order

// Make a wave to contain an exponential with decay constant pw[3]. The 
// wave needs enough points to allow any reasonable decay constant to 
// get really close to zero. The scaling is made symmetric about zero to 
// avoid an X offset from Convolve/A

// Make power decay wave

// add at least 15 lifetimes onto end of wave to deal with end effects 
// so we can delete them later
Variable endp = dimsize(xw,0)
Variable xfinal = xw[(endp-1)] + 15*(1/pw[3])

Variable PFdeta = (xfinal-xw[0])/endp
Variable NdelPts = round(15*(1/pw[3])/PFdeta)

//Variable nYPnts = 10*max(resolutionFactor*numpnts(yw), nexpWavePnts) 
//Variable nYPnts = 1*max(numpnts(yw), nexpWavePnts) 
Variable nYPnts = numpnts(yw)
Make/D/O/N=(nYPnts) powerwave

setscale/I x xw[0],xfinal, powerwave
Variable dtpw
dtpw = (xfinal - xw[0])/nYPnts
powerwave = ((x-pw[1])>0)*(((x-pw[1])^(-pw[4])))/(dtpw^(-pw[4])) 

//Variable dT = deltax(yw)
// Arbitrary choice, make sure we are oversmapling exp rise
//Variable resolutionFactor = 100
//Variable dT =(1/(resolutionFactor*pw[3]))

// Calculate suitable number points for the gaussian. Length of 
//  wave is doubled so function can start 
// in the middle; +1 to make it odd so  starts at t=0, and t=0 is 
// exactly the middle point. That is better for the convolution.

// make sure there are at least 15 lifetimes being represented in the wave
// this gets us to e^-15 = 3 *10^-7
//Variable nexpWavePnts = round((15*(1/pw[3]))/dT)*2 
// make sure wave scaling is identical to above
Variable nexpWavePnts =round((15*(1/pw[3]))/dtpw)*2
Make/D/O/N=(nexpWavePnts) expwave	// long enough to allow decay to zero 
Setscale/P x -dtpw*(nexpWavePnts/2),dtpw,expwave

// fill expwave with exponential decay 
//expwave = ((x-pw[1])>=0)*(dT/pw[3])*exp(-pw[3]*(x-pw[1]))
expwave = ((x-pw[1])>=0)*exp(-pw[3]*(x-pw[1]))

// Normalize exponential so that convolution doesn't change 
// the amplitude of the result 
//Variable sumexp 
//sumexp = sum(expwave, -inf,inf)
//expwave /= sumexp

// remove NaNs
powerwave = NumType(powerwave) == 2 ? 0:powerwave

// duplicate before convolution for trouble shooting
Duplicate/O powerwave, PWo

// and convolve with the exponential; NOTE /A 
convolve/A expwave, powerwave

// remove end effects
Variable StDel = DimSize(powerwave,0)-NdelPts
DeletePoints StDel, NdelPts, powerwave
yw = pw[0]+pw[2]*powerwave(xw[p])

End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function ExpConvPwEG(pw, yw, xw) : FitFunc 
WAVE pw, yw, xw
// pw[0] = baseline offset
// pw[1] = timezero
// pw[2] = amplitude 
// pw[3] = exponential rate constant
// pw[4] = power order
// pw[5] = ExpGauss Ampl
// pw[6] = pulsewidth
// pw[7] = expgauss rate

// Make a wave to contain an exponential with decay constant pw[3]. The 
// wave needs enough points to allow any reasonable decay constant to 
// get really close to zero. The scaling is made symmetric about zero to 
// avoid an X offset from Convolve/A

// determine shortest necessary timescale
Variable Gmax
Gmax = max(pw[3], pw[7])

Variable dtpw

// Make power decay wave

// add at least 15 lifetimes onto end of wave to deal with end effects 
// so we can delete them later
Variable endp = dimsize(xw,0)
Variable xfinal = xw[(endp-1)] + 15*(1/pw[3])

Variable PFdeta = (xfinal-xw[0])/endp
Variable NdelPts = round(15*(1/pw[3])/PFdeta)

//Variable nYPnts = 10*max(resolutionFactor*numpnts(yw), nexpWavePnts) 
//Variable nYPnts = 1*max(numpnts(yw), nexpWavePnts) 
Variable nYPnts = numpnts(yw)
// make sure dT oversamples this rate by a factor of 10 below
dtpw = min((xfinal - xw[0])/nYPnts, Gmax/10)
// redefine nYPnts
nYPnts = (xfinal - xw[0])/dtpw

Make/D/O/N=(nYPnts) powerwave
Make/D/O/N=(nYPnts) EGwave

setscale/I x xw[0],xfinal, powerwave
setscale/I x xw[0],xfinal, EGwave

dtpw = (xfinal - xw[0])/nYPnts
powerwave = ((x-pw[1])>0)*(((x-pw[1])^(-pw[4])))/(dtpw^(-pw[4])) 
EGwave = pw[5]*ExpGauss((x-pw[1]),pw[7],pw[6])

//Variable dT = deltax(yw)
// Arbitrary choice, make sure we are oversmapling exp rise
//Variable resolutionFactor = 100
//Variable dT =(1/(resolutionFactor*pw[3]))

// Calculate suitable number points for the gaussian. Length of 
//  wave is doubled so function can start 
// in the middle; +1 to make it odd so  starts at t=0, and t=0 is 
// exactly the middle point. That is better for the convolution.

// make sure there are at least 15 lifetimes being represented in the wave
// this gets us to e^-15 = 3 *10^-7
//Variable nexpWavePnts = round((15*(1/pw[3]))/dT)*2 
// make sure wave scaling is identical to above
Variable nexpWavePnts =round((15*(1/pw[3]))/dtpw)*2
Make/D/O/N=(nexpWavePnts) expwave	// long enough to allow decay to zero 
Setscale/P x -dtpw*(nexpWavePnts/2),dtpw,expwave

// fill expwave with exponential decay 
//expwave = ((x-pw[1])>=0)*(dT/pw[3])*exp(-pw[3]*(x-pw[1]))
expwave = ((x-pw[1])>=0)*exp(-pw[3]*(x-pw[1]))

// Normalize exponential so that convolution doesn't change 
// the amplitude of the result 
//Variable sumexp 
//sumexp = sum(expwave, -inf,inf)
//expwave /= sumexp

// remove NaNs
powerwave = NumType(powerwave) == 2 ? 0:powerwave

// duplicate before convolution for trouble shooting
Duplicate/O powerwave, PWo

// and convolve with the exponential; NOTE /A 
convolve/A expwave, powerwave

// add expgauss data
powerwave = pw[0]+pw[2]*powerwave+EGwave

// remove end effects
Variable StDel = DimSize(powerwave,0)-NdelPts
DeletePoints StDel, NdelPts, powerwave
yw = powerwave(xw[p])

End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function ExpConvPwEG2(pw, yw, xw) : FitFunc 
WAVE pw, yw, xw
// pw[0] = baseline offset
// pw[1] = timezero
// pw[2] = amplitude 
// pw[3] = exponential rate constant
// pw[4] = power order
// pw[5] = ExpGauss Ampl 1
// pw[6] = pulsewidth
// pw[7] = expgauss rate 1
// pw[8] = EG amp 2
// pw[9] = EG rate 2

// Make a wave to contain an exponential with decay constant pw[3]. The 
// wave needs enough points to allow any reasonable decay constant to 
// get really close to zero. The scaling is made symmetric about zero to 
// avoid an X offset from Convolve/A

// determine shortest necessary timescale
Variable Gmax
Gmax = max(pw[8], pw[9])
Gmax = max(pw[3], pw[9])

Variable dtpw

Variable endp = dimsize(xw,0)
// add at least 15 lifetimes onto end of wave to deal with end effects 
// so we can delete them later
Variable xfinal = xw[(endp-1)] + 15*(1/pw[3])
Variable PFdeta = (xfinal-xw[0])/endp
Variable NdelPts = round(15*(1/pw[3])/PFdeta)

//Variable nYPnts = 10*max(resolutionFactor*numpnts(yw), nexpWavePnts) 
//Variable nYPnts = 1*max(numpnts(yw), nexpWavePnts) 
Variable nYPnts = numpnts(yw)
// make sure dT oversamples this rate by a factor of 10 below
dtpw = min((xfinal - xw[0])/nYPnts, Gmax/10)
// redefine nYPnts
nYPnts = (xfinal - xw[0])/dtpw
Make/D/O/N=(nYPnts) powerwave
Make/D/O/N=(nYPnts) EGwave

// Make power decay wave

setscale/I x xw[0],xfinal, powerwave
setscale/I x xw[0],xfinal, EGwave

powerwave = ((x-pw[1])>0)*(((x-pw[1])^(-pw[4])))/(dtpw^(-pw[4])) 
EGwave = pw[5]*ExpGauss((x-pw[1]),pw[7],pw[6]) + pw[8]*ExpGauss((x-pw[1]),pw[9],pw[6])

//Variable dT = deltax(yw)
// Arbitrary choice, make sure we are oversmapling exp rise
//Variable resolutionFactor = 100
//Variable dT =(1/(resolutionFactor*pw[3]))

// Calculate suitable number points for the gaussian. Length of 
//  wave is doubled so function can start 
// in the middle; +1 to make it odd so  starts at t=0, and t=0 is 
// exactly the middle point. That is better for the convolution.

// make sure there are at least 15 lifetimes being represented in the wave
// this gets us to e^-15 = 3 *10^-7
//Variable nexpWavePnts = round((15*(1/pw[3]))/dT)*2 
// make sure wave scaling is identical to above
Variable nexpWavePnts =round((15*(1/pw[3]))/dtpw)*2
Make/D/O/N=(nexpWavePnts) expwave	// long enough to allow decay to zero 
Setscale/P x -dtpw*(nexpWavePnts/2),dtpw,expwave

// fill expwave with exponential decay 
//expwave = ((x-pw[1])>=0)*(dT/pw[3])*exp(-pw[3]*(x-pw[1]))
expwave = ((x-pw[1])>=0)*exp(-pw[3]*(x-pw[1]))

// Normalize exponential so that convolution doesn't change 
// the amplitude of the result 
//Variable sumexp 
//sumexp = sum(expwave, -inf,inf)
//expwave /= sumexp

// remove NaNs
powerwave = NumType(powerwave) == 2 ? 0:powerwave

// duplicate before convolution for trouble shooting
Duplicate/O powerwave, PWo

// and convolve with the exponential; NOTE /A 
convolve/A expwave, powerwave

// add expgauss data
powerwave = pw[0]+pw[2]*powerwave+EGwave

// remove end effects
Variable StDel = DimSize(powerwave,0)-NdelPts
DeletePoints StDel, NdelPts, powerwave
yw = powerwave(xw[p])

End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function ExpConvPwEG2c(pw, yw, xw) : FitFunc 
WAVE pw, yw, xw
// pw[0] = baseline offset
// pw[1] = timezero
// pw[2] = amplitude 
// pw[3] = exponential rate constant
// pw[4] = power order
// pw[5] = ExpGauss Ampl 1
// pw[6] = pulsewidth
// pw[7] = EG amp 2
// pw[8] = EG rate 2

// EG rate 1 fixed to pw[3]

// Make a wave to contain an exponential with decay constant pw[3]. The 
// wave needs enough points to allow any reasonable decay constant to 
// get really close to zero. The scaling is made symmetric about zero to 
// avoid an X offset from Convolve/A

// determine shortest necessary timescale
Variable Gmax
//Gmax = max(pw[8], pw[9])
Gmax = max(pw[3], pw[8])

Variable dtpw

Variable endp = dimsize(xw,0)
// add at least 15 lifetimes onto end of wave to deal with end effects 
// so we can delete them later
Variable xfinal = xw[(endp-1)] + 15*(1/pw[3])
Variable PFdeta = (xfinal-xw[0])/endp
Variable NdelPts = round(15*(1/pw[3])/PFdeta)

//Variable nYPnts = 10*max(resolutionFactor*numpnts(yw), nexpWavePnts) 
//Variable nYPnts = 1*max(numpnts(yw), nexpWavePnts) 
Variable nYPnts = numpnts(yw)
// make sure dT oversamples this rate by a factor of 10 below
dtpw = min((xfinal - xw[0])/nYPnts, Gmax/10)
// redefine nYPnts
nYPnts = (xfinal - xw[0])/dtpw
Make/D/O/N=(nYPnts) powerwave
Make/D/O/N=(nYPnts) EGwave

// Make power decay wave

setscale/I x xw[0],xfinal, powerwave
setscale/I x xw[0],xfinal, EGwave

powerwave = ((x-pw[1])>0)*(((x-pw[1])^(-pw[4])))/(dtpw^(-pw[4])) 
EGwave = pw[5]*ExpGauss((x-pw[1]),pw[3],pw[6]) + pw[7]*ExpGauss((x-pw[1]),pw[8],pw[6])

//Variable dT = deltax(yw)
// Arbitrary choice, make sure we are oversmapling exp rise
//Variable resolutionFactor = 100
//Variable dT =(1/(resolutionFactor*pw[3]))

// Calculate suitable number points for the gaussian. Length of 
//  wave is doubled so function can start 
// in the middle; +1 to make it odd so  starts at t=0, and t=0 is 
// exactly the middle point. That is better for the convolution.

// make sure there are at least 15 lifetimes being represented in the wave
// this gets us to e^-15 = 3 *10^-7
//Variable nexpWavePnts = round((15*(1/pw[3]))/dT)*2 
// make sure wave scaling is identical to above
Variable nexpWavePnts =round((15*(1/pw[3]))/dtpw)*2
Make/D/O/N=(nexpWavePnts) expwave	// long enough to allow decay to zero 
Setscale/P x -dtpw*(nexpWavePnts/2),dtpw,expwave

// fill expwave with exponential decay 
//expwave = ((x-pw[1])>=0)*(dT/pw[3])*exp(-pw[3]*(x-pw[1]))
expwave = ((x-pw[1])>=0)*exp(-pw[3]*(x-pw[1]))

// Normalize exponential so that convolution doesn't change 
// the amplitude of the result 
//Variable sumexp 
//sumexp = sum(expwave, -inf,inf)
//expwave /= sumexp

// remove NaNs
powerwave = NumType(powerwave) == 2 ? 0:powerwave

// duplicate before convolution for trouble shooting
Duplicate/O powerwave, PWo

// and convolve with the exponential; NOTE /A 
convolve/A expwave, powerwave

// add expgauss data
powerwave = pw[0]+pw[2]*powerwave+EGwave

// remove end effects
Variable StDel = DimSize(powerwave,0)-NdelPts
DeletePoints StDel, NdelPts, powerwave
yw = powerwave(xw[p])

End

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////


Function ExpConvPwEG3c(pw, yw, xw) : FitFunc 
WAVE pw, yw, xw
// pw[0] = baseline offset
// pw[1] = timezero
// pw[2] = amplitude 
// pw[3] = exponential rate constant
// pw[4] = power order
// pw[5] = ExpGauss Ampl 1
// pw[6] = pulsewidth
// pw[7] = EG amp 2
// pw[8] = EG rate 2
// pw[9] = EG amp 3
// pw[10] = EG rate 3

// EG rate 1 fixed to pw[3]

// Make a wave to contain an exponential with decay constant pw[3]. The 
// wave needs enough points to allow any reasonable decay constant to 
// get really close to zero. The scaling is made symmetric about zero to 
// avoid an X offset from Convolve/A

// determine shortest necessary timescale
Variable Gmax
//Gmax = max(pw[8], pw[9])
Gmax = max(max(pw[3], pw[8]), pw[10])

Variable dtpw

Variable endp = dimsize(xw,0)
// add at least 15 lifetimes onto end of wave to deal with end effects 
// so we can delete them later
Variable xfinal = xw[(endp-1)] + 15*(1/pw[3])
Variable PFdeta = (xfinal-xw[0])/endp
Variable NdelPts = round(15*(1/pw[3])/PFdeta)

//Variable nYPnts = 10*max(resolutionFactor*numpnts(yw), nexpWavePnts) 
//Variable nYPnts = 1*max(numpnts(yw), nexpWavePnts) 
Variable nYPnts = numpnts(yw)
// make sure dT oversamples this rate by a factor of 10 below
dtpw = min((xfinal - xw[0])/nYPnts, Gmax/10)
// redefine nYPnts
nYPnts = (xfinal - xw[0])/dtpw
Make/D/O/N=(nYPnts) powerwave
Make/D/O/N=(nYPnts) EGwave

// Make power decay wave

setscale/I x xw[0],xfinal, powerwave
setscale/I x xw[0],xfinal, EGwave

powerwave = ((x-pw[1])>0)*(((x-pw[1])^(-pw[4])))/(dtpw^(-pw[4])) 
EGwave = pw[5]*ExpGauss((x-pw[1]),pw[3],pw[6]) + pw[7]*ExpGauss((x-pw[1]),pw[8],pw[6])+ pw[9]*ExpGauss((x-pw[1]),pw[10],pw[6])

//Variable dT = deltax(yw)
// Arbitrary choice, make sure we are oversmapling exp rise
//Variable resolutionFactor = 100
//Variable dT =(1/(resolutionFactor*pw[3]))

// Calculate suitable number points for the gaussian. Length of 
//  wave is doubled so function can start 
// in the middle; +1 to make it odd so  starts at t=0, and t=0 is 
// exactly the middle point. That is better for the convolution.

// make sure there are at least 15 lifetimes being represented in the wave
// this gets us to e^-15 = 3 *10^-7
//Variable nexpWavePnts = round((15*(1/pw[3]))/dT)*2 
// make sure wave scaling is identical to above
Variable nexpWavePnts =round((15*(1/pw[3]))/dtpw)*2
Make/D/O/N=(nexpWavePnts) expwave	// long enough to allow decay to zero 
Setscale/P x -dtpw*(nexpWavePnts/2),dtpw,expwave

// fill expwave with exponential decay 
//expwave = ((x-pw[1])>=0)*(dT/pw[3])*exp(-pw[3]*(x-pw[1]))
expwave = ((x-pw[1])>=0)*exp(-pw[3]*(x-pw[1]))

// Normalize exponential so that convolution doesn't change 
// the amplitude of the result 
//Variable sumexp 
//sumexp = sum(expwave, -inf,inf)
//expwave /= sumexp

// remove NaNs
powerwave = NumType(powerwave) == 2 ? 0:powerwave

// duplicate before convolution for trouble shooting
Duplicate/O powerwave, PWo

// and convolve with the exponential; NOTE /A 
convolve/A expwave, powerwave

// add expgauss data
powerwave = pw[0]+pw[2]*powerwave+EGwave

// remove end effects
Variable StDel = DimSize(powerwave,0)-NdelPts
DeletePoints StDel, NdelPts, powerwave
yw = powerwave(xw[p])

End