Function master_kin_fit(pw, yw, xw) : FitFunc
	Wave pw, yw, xw
	
	// pw[0] = temp (free parameter for empirical fit)
	// pw[1] = mass
	// pw[2] = kin_num (hold =1, =2, or =3)
	// pw[3] = branch_par (hold =0.5 for no splitting or kin_num 1)
	// pw[4] = A
	// pw[5] = tau1
	// pw[6] = tau2a (hold any for kin_num 1)
	// pw[7] = tau2b (hold any for no splitting or kin_num 1)
	// pw[8] = x0 (laser fire time: free parameter for empirical fit)
	// pw[9] = y0
	// pw[10] = u (speed parameter, for supersonic fit. hold =0 unless otherwise noted)
	
	variable term1, term2, term3, term4, term5
	variable normarea1, normarea2, stepsize = xw[1]-xw[0]
	variable resolution = 1 // THIS ISN'T ARBITRARY!! it's the t bin width
	
	Make/Free/D/N=(numpnts(yw)*resolution) effusion_filter, kin_func, heaviside
	SetScale/P x 0,stepsize/resolution,"", effusion_filter // not generalised yet, sorry!
	SetScale/P x xw[0],stepsize/resolution,"", kin_func, heaviside
	
	term1 = (-pw[1]*1.66053906660e-4*0.025^2/(2*1.380649*pw[0]))
	effusion_filter = exp(term1*((1/(x*0.001)^2)-pw[10]/(0.025^2)))/(x*0.001)^4
	effusion_filter = NumType(effusion_filter)==2 ? 0 : effusion_filter // Replace NaNs with zero
	normarea1 = area(effusion_filter)
	effusion_filter = effusion_filter/normarea1	// normalise area to 1
	
	if(pw[2] == 1)								// single exponential with x, y offset
		kin_func = pw[4]*exp(-(x-pw[8])/pw[5])+pw[9]
	elseif(pw[2] == 2)							// triple exponential with zero start, then x, y offset
		term1 = ((1/2)+pw[3])*(pw[6]/(pw[5]-pw[6]))
		term2 = ((1/2)-pw[3])*(pw[7]/(pw[5]-pw[7]))
		term3 = term1+term2
		kin_func = (term3*exp(-(x-pw[8])/pw[5])-term1*exp(-(x-pw[8])/pw[6])-term2*exp(-(x-pw[8])/pw[7]))*pw[4]+pw[9]
	elseif(pw[2] == 3)							// triple exponential with zero derivative at start, with x, y offset
		term1 = ((1/2)+pw[3])*(pw[6]/(pw[5]-pw[6]))
		term2 = ((1/2)-pw[3])*(pw[7]/(pw[5]-pw[7]))
		term3 = ((1/2)+pw[3])*(pw[5]/(pw[5]-pw[6]))
		term4 = ((1/2)-pw[3])*(pw[5]/(pw[5]-pw[7]))
		term5 = term3+term4
		kin_func = (1-term5*exp(-(x-pw[8])/pw[5])+term1*exp(-(x-pw[8])/pw[6])+term2*exp(-(x-pw[8])/pw[7]))*pw[4]+pw[9]	
	else
		print("this needs to be 1, 2 or 3")
	endif
	
	heaviside = ((1+abs(x-pw[8])/(x-pw[8]))/2)	
	heaviside = NumType(heaviside)==2 ? 0.5 : heaviside // Replace NaNs with one half for a heaviside
	kin_func  = NumType(kin_func)==2  ? 0   : kin_func  // Replace NaNs with zero
	
	// the heaviside is to zero before the laser fires, since the function is nonzero outside the boundary conditions otherwise
	kin_func = kin_func*heaviside
	normarea1 = area(kin_func)
	convolve effusion_filter kin_func
	normarea2 = area(kin_func)
	kin_func = kin_func*normarea1/normarea2

	if (resolution!=1)
		yw = mean(kin_func,pnt2x(kin_func,p*resolution),pnt2x(kin_func,(p+1)*resolution))
	else
		yw = kin_func
	endif
End