function [a,aerr,chisq,yfit,corr] = levmar(x,y,sig,fitfun,a0,dYda,sgn)

% levmar.m Fits a nonlinear function to data using the Marquardt
%                     method discussed in Bevington and Robinson in Section 8.6.
%    [a,aerr,chisq,yfit] = levmar(x,y,sig,fitfun,a0,dYda,sgn)
%
%    Inputs:  x -- the x data to fit
%             y -- the y data to fit
%             sig -- the uncertainties on the data points
%             fitfun -- the name of the function to fit to
%             a0 -- the initial guess at the parameters
%             sgn, dYda -- if sgn=0, use numerical derivatives;
%                          if sgn=1, use analytic expressions for
%                          derivatives stored in the cell "dYda"
%    Outputs: a -- the best fit parameters
%             aerr -- the errors on these parameters
%             chisq -- the final value of chi-squared
%             yfit -- the value of the fitted function
%                     at the points in x
%             corr -- the error matrix
%*****************************************
%*** Parameters you may need to modify ***
%*****************************************
stepsize  = abs(a0)*0.001;    % the step of increase of a when calculating the derivatives
chicut = 0.01;                % Maximum differential allowed between successive chi^2 values
silent = 0;                   % Change to 1 to supress output.

%*****************************************

a = a0;
chi2 = calcchi2(x,y,sig,fitfun,a);     % computer chi square as in step 1
chi1 = chi2+chicut*2;
lambda = 0.001;                        % start lambda at 0.001
[dum,nparm] = size(a);

i=0;
while (abs(chi2-chi1))>chicut
   i=i+1;
  if silent~=1
   fprintf(1,'a = ');
   fprintf(1,'%6.2f ',a);
   fprintf(1,'\t ChiSq = %6.2f', chi2);
   fprintf(1,'\t diff = %6.2f\n', abs(chi2-chi1));
   fprintf(1,'\n');
  end

 chinew = chi2 + 1;
 while chinew>chi2+chicut
     deltaa = calcdeltaa(x,y,sig,fitfun,a,stepsize,lambda,dYda,sgn);    
     anew   = a + deltaa;
     chinew = calcchi2(x,y,sig,fitfun,anew);
     if chinew>chi2
         lambda = lambda*10;                     
     end
 end

 lambda = lambda/10;   
 a      = anew;         
 chi1   = chi2;
 chi2   = chinew;
end

if silent~=1
   fprintf(1,'a = ');
   fprintf(1,'%6.2f ',a);
   fprintf(1,'\t ChiSq = %6.2f', chi2);
   fprintf(1,'\t diff = %6.2f\n', abs(chi2-chi1));
end
i

% calculate the returned values
% the diagonal terms of the error matrix are the variances 
corr = calcinvalpha(x,y,sig,fitfun,a,stepsize,lambda,dYda,sgn)
for i=1:nparm
   aerr(i) = sqrt(corr(i,i));
end
chisq = calcchi2(x,y,sig,fitfun,a);
yfit = feval(fitfun,x,a);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this function just calculates the value of chi^2
function chi2 = calcchi2(x,y,sig,fitfun,a)

chi2 = sum( ((y-feval(fitfun,x,a)) ./sig).^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this function calculates the value of deltaa according to the current
% choice of a. See bevington p. 157 equation (8.28)
function deltaa = calcdeltaa(x,y,sig,fitfun,a,stepsize,lambda,dYda,sgn)

[dum,nparm] = size(a);
[ndata,dum] = size(x);
beta = zeros(1,nparm);

corr = calcinvalpha(x,y,sig,fitfun,a,stepsize,lambda,dYda,sgn);
der = calcderiv(x,y,sig,fitfun,a,stepsize,dYda,sgn);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this loop calculates the vector beta see Bevington Page 293, Program8.6
for k=1:nparm
   for i=1:ndata
       beta(k) = beta(k)+(y(i)-feval(fitfun,x(i),a))*der(i,k)/sig(i)/sig(i);
   end
end
deltaa = beta*corr;	% Bevington Equation 8.28

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this function calculates the square matrix alpha and its inverse corr
% See Bevington p. 159 eq.(8.35) and Page 293, Program 8.6

function corr = calcinvalpha(x,y,sig,fitfun,a,stepsize,lambda,dYda,sgn);

[dum,nparm] = size(a);
[ndata,dum] = size(x);
alpha = zeros(nparm,nparm);

der = calcderiv(x,y,sig,fitfun,a,stepsize,dYda,sgn);
for j=1:nparm
   for k=1:nparm
       for i=1:ndata
           alpha(j,k) = alpha(j,k)+der(i,j)*der(i,k)/sig(i)/sig(i);
       end
   end
end
for j=1:nparm    alpha(j,j) = (1+lambda)*alpha(j,j);
end
corr = inv(alpha);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%this function calculates the values of dY[i]/da[j] for all i and j
function der = calcderiv(x,y,sig,fitfun,a,stepsize,dYda,sgn)

[dum,nparm] = size(a);
[ndata,dum] = size(x);
der = zeros(ndata,nparm);

% if sgn=0, use numerival derivatives
if sgn==0
   for i=1:ndata
       y0 = feval(fitfun,x(i),a);
       for j=1:nparm
           a(j) = a(j) + stepsize(j);
           y1 = feval(fitfun,x(i),a);
           a(j) = a(j) - stepsize(j);
           der(i,j) = (y1-y0)/stepsize(j);	% Bevington Equation A.24
       end
   end
else   % if sgn=1, use the analytic expressions for derivatives
   for i=1:ndata
       for j=1:nparm
           der(i,j) = feval(dYda{j},x(i),a);
       end
   end
end
