function [pars,parscovar]=fitLorentz(x,yexp,yerrors,guess,minBound,maxBound)
% Fits the experimental data by a Lorentz function
% 
% Václav Římal, 9. 12. 2019
% 
% Lorentz is described by parameters [I x0 g], where I is the intensity,
% x0 is the median and g is the FWHM
% 
% yerrors: optional
% guess: initial set of parameters
% minBound, maxBound: optional, specify the lower and upper limits of
% parameters
% 
% returns the optimal parameters and their covariance matrix
% 
% Example usage:
% [p,c]=fitLorentz(x,y,[],[5 0 0.5])
% 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% are the yerrors included?
% plot the initial state
if (nargin<3)||isempty(yerrors)
	plot(x,yexp,'x');
    yerrors=ones(size(yexp));
else
	errorbar(x,yexp,yerrors,'x');
end
hold on;
plot(x,Lorentz(x,guess),'--');
xlabel('\delta (^{39}K) / ppm');


% weighting of data by their errors is made by
% substitution yexpNove = 0 and yfitNove=(yfit-yexp)./yerrors
fce=@(p,x)weighteddifference(yexp,Lorentz(x,p),yerrors);


% set some options for nonlin fitting procedure
options=optimset('Display','final',...
        'MaxFunEvals',5000,'MaxIter',300,'TolFun',1e-11,'TolX',1e-12);       


% fitting - minimization of chi-square by lsqcurvefit
if (nargin<5)||isempty(minBound)||(nargin<6)||isempty(maxBound)
        [pars,fitinfo,res,exflag,output,lambda,jacobian]=lsqcurvefit(...
            fce,guess,x,zeros(size(yexp)),[],[],options);
else
        [pars,fitinfo,res,exflag,output,lambda,jacobian]=lsqcurvefit(...
            fce,guess,x,zeros(size(yexp)),minBound,maxBound,options);
end    
% calculate errors from Jacobian
% weighting by experimental errors is done already by the
% weighteddifference function used as the objective function by
% lsqcurvefit
       
parscovar=full((jacobian'*jacobian)^(-1));
        

% units are as follows: 
% [jacobian] = [y]/[yerrors]/[pars] = 1/[pars];
% [parscovar] = [jacobian]^(-2) = [pars]^2;

        
% if real yerrors unavailable, then:
if isempty(yerrors) || ~any(isfinite(yerrors))
	% covariance weighting by residuals    
	% normres = norm(res)/sqrt(df)    
	normres=norm(res)/sqrt( numel(yexp)-size(pars,1));    
	parscovar = parscovar*normres^2;    
end
         
plot(x,Lorentz(x,pars));
hold off;

legend({'experiment';'initial guess';'best fit'});
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definition of the Lorentz function
% 
function y=Lorentz(x,p)
	I=p(1);
	x0=p(2);
	g=p(3);

	y=I*g./(g*g+(x-x0).^2)/pi;
end

