function [pout, yfit, res] = fitter(x,y,f, pin)
	% x - x data points vector
	% y - y data points vector
	% f - handle to function f(x,p) to be fitted, i.e. model
	% pin - function/model parameters initial guess vector 
	%
	% pout - parameters vector providing the best fit
	% yfit - y values calculated by the functions f(x,pout)
	% res - residuals (y-yfit)

	function yfit=yfitted(x,p)
		f_of_x=@(x) f(x,p);
		yfit=arrayfun(f_of_x,x);
	end

	function res=residuals(x,p)
		yfit = yfitted(x,p);
		res = y-yfit;
	end
		
	% our merit function sum of residuals squared
	function E=chi_sq(p)
		res=residuals(x,p);
		E = sum( (res).^2 );
	end

	% the main job of fitting 
	% (i.e.  optimization/best parameter finding)
    % is done by fminsearch
	pout=fminsearch(@chi_sq, pin); %fit is complete

	yfit=yfitted(x,pout);
	res=residuals(x,pout);
end

