function f = ezfit(varargin)
%EZFIT   Fit data with arbitrary fitting function
%   EZFIT(FUN) fits the active curve with the function FUN. See below for
%   the syntax of FUN. If FUN is not specified, 'linear' is used.
%
%   By default, the first curve in the active figure is fitted - see
%   FITPARAM to change this default behavior. To fit another curve, select
%   it before calling EZFIT.  If some data are selected by the "Data
%   Brushing" tool (only for Matlab >= 7.6), only those data are fitted.
%
%   EZFIT(X,Y,FUN) or EZFIT(Y,FUN) fits the data (X,Y) (or Y) using the
%   function FUN. X and Y must be vectors of equal length. If X is not
%   specified, X=[1, 2, 3...] is assumed.
%
%   EZFIT(X,Y,FUN), where X is a 1-by-N vector and Y is a 2-by-N matrix,
%   also specifies the weights for Y(2,:) (error bars). By default, when
%   Y is a 1-by-N vector, all the weights are 1.
%
%   Note that EZFIT only computes the coefficients, but does not display the
%   fit. Use SHOWFIT to display the fit.
%
%   The function string FUN can be:
%      - the name of a predefined fitting function (see below).
%      - the name of a user-defined fitting function (see EDITFIT).
%      - an equation, in the form 'y(x)=...', where 'x' represents the
%        X-data, and all the other variables are parameters to be fitted
%        ('a', 'x_0', 'tau', ...). Example: 'y(x)=a*sin(b*x)'. If the
%        left-hand-side 'y(x)' is not specified, 'x' is taken for the
%        X-Data. All the parameter names are accepted, except Matlab
%        reserved strings ('sin', 'pi', ...)
%
%   The predefined fitting functions are:
%      - linear             y = m * x
%      - affine or poly1    y = a*x + b
%      - poly{n}            y = a0 + a1 * x + ... + an * x^n
%      - power              y = c*x^n
%      - sin                y = a * sin (b * x)
%      - cos                y = a * cos (b * x)
%      - exp                y = a * exp (b * x)
%      - log                y = a * log (b * x)
%      - cngauss            y = exp(-x^2/(2*s^2))/(2*pi*s^2)^(1/2)
%      - cfgauss            y = a*exp(-x^2/(2*s^2))
%      - ngauss             y = exp(-(x-x0)^2/(2*s^2))/(2*pi*s^2)^(1/2)
%      - gauss              y = a*exp(-(x-x0)^2/(2*s^2))
%   'ngauss' is a 2-parameters normalized Gaussian, and 'gauss' is a
%   3-parameters non-normalized (free) Gaussian. 'cngauss' and 'cfgauss'
%   are centered normalized and centered free Gaussian, respectively.
%
%   EZFIT is based on Matlab's built-in FMINSEARCH function (Nelder-Mead
%   method), which performs an unconstrained nonlinear minimization of
%   the SSR (sum of squared residuals) with respect to the various
%   parameters. 
%
%   The correlation coefficient R is defined as (SSR/(SSE+SSR))^(1/2), where
%      SSR = sum ((y_fit - mean(y)).^2)   % sum of squared residuals
%      SSE = sum ((y_fit - y).^2)         % sum of squared errors
%   (see http://mathworld.wolfram.com/CorrelationCoefficient.html)
%
%   Nonlinear minimization requires starting guesses (or starting estimates)
%   for the fit parameters. By default, all the starting guesses are taken
%   as 1, or, when using predefined fits (e.g., exp, gauss, power...), the
%   starting guesses are determined depending on the range of the data to
%   be fitted. However, in most cases, values closer to the expected
%   result should be specified to "help" the convergence. It is sufficient
%   to choose values that have the correct sign and correct order of
%   magnitude, e.g. 0.01, 1, 100...
%
%   The starting guesses for the parameters may be specified in two ways:
%     - directly in the string FUN, after the fit definition:
%          'c0 + a*sin(pi*x/lambda); c0=1; a=0.1; lambda=100'
%          ('!' or '$' may also be used instead of ';').
%     - by specifying them as an additional input argument for EZFIT:
%          EZFIT(x,y,'c0 + a*sin(pi*x/lambda)',[0.1 1 100]);
%       Note that in this case the parameters must be ordered alphabetically.
%   Note that if both methods are used, only the starting guesses given in
%   the string FUN are considered.
%
%   By default, Y is fitted in linear mode. If you want to fit LOG(Y)
%   instead, you must specify the option 'log' to the string FUN, separated
%   by the symbol ';' or '$' or '!' (eg, FUN='a*x^n;log'). This is
%   specially useful to fit power laws with equally weighted points in a
%   log scale.  If nothing specified, the option 'lin' is used.
%
%   Example:
%        plotsample('power')
%   and compare the results of:
%        ezfit('power;lin')
%        ezfit('power;log')
%
%   F = EZFIT(...) also returns a structure F having the following fields:
%      - name       name of the fit
%      - eq         equation of the fit
%      - param      cell array of strings: names of the parameters
%      - m          values of the coefficients
%      - m0         initial guess for the coefficients
%      - r          correlation coefficient R (Pearson's correlation)
%      - fitmode    'lin' (y is fitted) or 'log' (log(y) is fitted) mode
%
%   This structure F can be further used with SHOWFIT, DISPEQFIT,
%   SHOWEQBOX, MAKEVARFIT and EDITCOEFF.
%
%   From F, you can get the values of the fitted parameters. If you want to
%   create in the current Matlab workspace the variables associated to
%   these parameters, use MAKEVARFIT (or set the option 'automakevarfit'
%   to 'on' in FITPARAM).
%
%   Examples:
%     plotsample damposc
%     f = ezfit('u(t) = c + u_a * sin(2*pi*t/T) * exp(-t/tau); T=5; tau=20');
%     showfit(f);
%     
%     plotsample poly2
%     [x,y] = pickdata;
%     f = ezfit(x, y, 'z(v) = poly3');
%     editcoeff(f);
%
%     plotsample poly2
%     f = ezfit('beta(z) = poly2');
%     showfit(f, 'fitcolor', 'red', 'fitlinewidth', 2);
%
%   See also SHOWFIT, PLOTSAMPLE, DISPEQFIT, EDITCOEFF,
%   FMINSEARCH, MAKEVARFIT, FITPARAM.


%   F. Moisy, moisy_at_fast.u-psud.fr
%   Revision: 2.51,  Date: 2011/10/04
%   This function is part of the EzyFit Toolbox


% History:
% 2005/05/12: v1.00, first version.
% 2005/05/20: v1.10, Use 'eval', with generic functions
% 2005/05/21: v1.11, added the 'lin','log' options.
% 2005/05/24: v1.12, option 'log' by default for 'power' and 'exp'.
% 2005/07/27: v1.13, cosmetics.
% 2005/09/03: v1.14, check arg.
% 2005/10/07: v1.15, gaussian fits added (ngauss and fgauss, centered/not)
% 2005/10/20: v1.16, help text changed.
% 2005/10/31: v1.20, also returns R. cste and poly{n} fits added. Initial
%                    guess defined within the fitting function string. The
%                    order of the output parameters is changed.
% 2005/11/05: v1.21, evaluate strings for initial guess in FUN
% 2005/12/06: v1.22, opens a dialog box if the polynomial order is not
%                    specified.
% 2006/01/13: v1.24, check for negative data in log mode
% 2006/01/19: v1.25, bug fixed from 1.24
% 2006/01/25: v1.26, check the matlab version
% 2006/02/08: v2.00, new syntax. The output argument is now a fit
%                    structure, and the fitting equation string accepts
%                    arbitrary parameter names.
% 2006/02/14: v2.10, lhs 'y(x)=...' now accepted. Now case sensitive
% 2006/03/07: v2.11, bug fixed from v1.25
% 2006/03/09: v2.20, weighted chi-square criterion (ie: error bars accepted
%                    for y) - undocumented
% 2006/09/04: v2.21, '!' and '$' may be used instead of ';' in the FUN
%                    string: this allows to pass the argument like:
%                    fit power!log
% 2006/09/28: v2.22, 'x^1' replaced by 'x' for 1st order polynomial
% 2006/10/18: v2.30, input raw and column vectors accepted. accepts
%                    additional input parameters to change the default settings.
% 2007/04/17: v2.31, help text improved; standard error messages
% 2007/05/16: v2.40, weigthed fits (v2.20) now documented, with bug fixed.
% 2007/08/18: v2.41, help text improved.
% 2007/09/17: v2.50, guess the initial guess m0 for predefined fits
% 2011/11/04: v2.51, bug fixed for fits in log coordinates
% 2012/06/28: v2.52, check version number removed (unstable)




% new v1.26, changed v2.11, removed v2.52
% if str2double(version('-release'))<14
%     error('Ezyfit:ezfit:compatibility','EzyFit requires Matlab 7 or higher.');
% end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fit parameters:  (new v2.30)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% loads the default fit parameters:
try
    fp=fitparam;
catch
    error('Ezyfit:ezfit:fitparamNotFound','''fitparam.m'' file not found.');
end


% change the default values of the fit parameters according to the
% additional input arguments:
for nopt=1:(nargin-1)
    if any(strcmp(varargin{nopt},fieldnames(fp))) % if the option string is one of the fit parameter
        fp.(varargin{nopt}) = varargin{nopt+1};
    end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input arguments:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==0, % if no input argument: do a linear fit of the current curve
    [x,y,h] = pickdata(fp);
    inputstr='linear';
elseif nargin==1
    if ischar(varargin{1}),    %    EZFIT('a*x+b')
        inputstr=varargin{1};
        [x,y,h] = pickdata(fp);
    else                       %    EZFIT(y)
        y=varargin{1};
        x=1:length(y);
        inputstr='linear';
    end
else   % 2 or more input arguments
    if ~isnumeric(varargin{1})           % EZFIT('fun',...)
        inputstr=varargin{1};
        [x,y,h] = pickdata(fp);
        if isnumeric(varargin{2})    % EZFIT('fun',m0,...)
            m0=varargin{2};
        end
    elseif (isnumeric(varargin{1}) && ~isnumeric(varargin{2}))       % EZFIT(y,'a*x+b',...)
        [y,inputstr]=deal(varargin{1:2});
        x=1:length(y);
        if nargin>2
            if isnumeric(varargin{3})   % EZFIT(y,'a*x+b',m0,...)    
                m0=varargin{3};       
            end
        end
    elseif (isnumeric(varargin{1}) && isnumeric(varargin{2}))   % EZFIT(x,y,...)
        [x,y]=deal(varargin{1:2});      
        if nargin>2
            if ischar(varargin{3})
                inputstr=varargin{3};      % EZFIT(x,y,'fun',...)
                if nargin>3
                    if isnumeric(varargin{4})      % EZFIT(x,y,'fun',m0,...)
                        m0=varargin{4};
                    end
                end
            else
                error('Ezyfit:ezfit:syntaxError','Syntax error. 3rd paramater of EZFIT should be a string.');
            end
        end
    end
end               



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some checks about x and y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% turn all input vectors into row (1xN) vectors (new v2.30):
if size(x,1)>size(x,2),
    x=x';
end
if size(y,1)>size(y,2),
    y=y';
end


% check for error bars for y (new v2.20, fixed 2.40)
if size(y,1)>1
    dy = y(2,:);  % second line = error bars (1/weight)
    y = y(1,:);   % first line = data
else
    dy = ones(1,length(y));    % default error bars: 1
end

if length(x)~=length(y),
    error('Ezyfit:ezfit:dimagree','X and Y dimensions must agree.');
end







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% processing of the string FUN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% cleans the input string:
inputstr = strrep(inputstr,' ','');
inputstr = strrep(inputstr,'!',';');
inputstr = strrep(inputstr,'$',';');
inputstr = strrep([inputstr ';'],';;',';'); % ensures that fun is terminated by a ';'.

% the name of the fit is by default the first part of the input string
p=findstr(inputstr,';'); p=p(1);
f.name = inputstr(1:(p-1));

% separates the first part (fitting function itself) and the remainder:
p = strfind(inputstr,';'); p=p(1);
fun = inputstr(1:(p-1));
remfun = inputstr((p+1):end);


% search for predefined fit or user-defined fit
usepredefinedfit = '';
[defaultfit, userfit] = loadfit;
for i=1:length(defaultfit),
    if strcmp(fun, defaultfit(i).name);
        fun = defaultfit(i).eq;
        usepredefinedfit = defaultfit(i).name;   % new v2.50
    end
end
for i=1:length(userfit),
    if strcmp(fun, userfit(i).name),
        fun = userfit(i).eq;
    end
end


% separates again the first part (fitting function itself) and the remainder:
% (because the predefined/user-defined part may itself contain ';')
fun = [strrep(fun,' ','') ';' remfun];
p=strfind(fun,';'); p=p(1);
remfun = fun((p+1):end);
fun = fun(1:(p-1));


% recognize if a lhs is present
peq = strfind(fun,'=');
if ~isempty(peq),
    lhs = fun(1:(peq-1));    % left-hand side
    rhs = fun((peq+1):end);  % right-hand side
else
    lhs = '';
    rhs = fun;
end


% process the lhs (if present)
if ~isempty(lhs),
    pob = strfind(lhs,'('); % position of opening bracket
    pcb = strfind(lhs,')'); % position of closing bracket
    if ~isempty(pob)
        if pob==1,
            f.yvar = 'y';
        else
            f.yvar = lhs(1:(pob-1));
        end
        f.xvar = lhs((pob+1):(pcb-1));
    else
        f.yvar = lhs;
        f.xvar = 'x';
    end
else % if no lhs present:
    f.xvar='x';
    f.yvar='y';
end


% process the 'poly' (polynomial fit) in the rhs
if strfind(rhs,'poly'),   % polynomial fit:
    order=str2double(rhs(5:end));
    if isempty(order), % new v1.22
        str_ord=inputdlg('Order of the polynomial fit','Polynomial order',1,{'2'});
        if ~isempty(str_ord),
            order=str2double(str_ord{1});
            f.name=['poly' str_ord{1}];
        else
            clear f;
            return;
        end
    end
    if order>20,
        error('Ezyfit:ezfit:invalidPolynomDegree','Invalid polynom degree.');
    end
    rhs = [fp.polynom_coeffname '0'];
    for i=1:order,
        if i>1
            rhs = [rhs '+' fp.polynom_coeffname num2str(i) '*' f.xvar '^' num2str(i)];
        else
            rhs = [rhs '+' fp.polynom_coeffname num2str(i) '*' f.xvar];   % new v2.22
        end
    end
end


% search for option 'lin' or 'log':
% (if several are present, use the last one)
% (if none is present, check the Y-scale of the current figure)
% (if no figure present, take 'lin').
fun = strrep([rhs ';' remfun], ';;', ';');
f.fitmode='';
plin=strfind(fun,';lin');
if ~isempty(plin), plin=plin(end); else plin=1; end
plog=strfind(fun,';log');
if ~isempty(plog), plog=plog(end); else plog=1; end
plast=max(plin,plog);
if plast==1, % if nothing specified
    f.fitmode='lin'; % use 'lin'
else
    f.fitmode=fun((plast+1):(plast+3));
end
fun=strrep(fun,';lin','');
fun=strrep(fun,';log','');


% check for negative data in log mode (new v1.23, fixed 1.24):
if (strcmp(f.fitmode,'log') && sum(y<=0)>0)
    disp('Warning: Zero or negative data ignored');
    nonzero=find(y>0);
    x=x(nonzero);
    y=y(nonzero);
end


% separates again the first part (fitting function itself) and the
% remainder:
p = strfind(fun,';'); p=p(1);
f.eq = fun(1:(p-1));
remfun = fun((p+1):end);


 % convert the equation in the matlab syntax (parameters named m(1),
 % m(2)...)
[eqml,param] = eq2ml(f.eq, f.xvar);
maxm = length(param); % number of parameters
if maxm==0  % new v2.20
    error('Ezyfit:ezfit:noParameter','No parameter to be fitted.');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% processing the initial guess
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% initial guess for the m(i) by default (all m(i)=1):
if ~exist('m0','var')
    m0=ones(1,maxm);
    switch usepredefinedfit   % new v2.50 : "magic" initial guesses
        case 'linear'   % m*x
            m0(1) = (y(end)-y(1))/(x(end)-x(1));
        case 'affine'   % a*x+b
            m0(1) = (y(end)-y(1))/(x(end)-x(1));  % a
            m0(2) = mean(y-m0(1)*x);              % b
        case 'affineshift'   % a*(x-b)
            m0(1) = (y(end)-y(1))/(x(end)-x(1));  % a
            m0(2) = mean(x-y/m0(1));              % b
        case 'power'    % a*x^n;
            % if both x and y are of constant sign:
            if sum(diff(sign(x)))==0 || sum(diff(sign(y)))==0
                m0(2) = log(y(end)/y(1)) / log(x(end)/x(1));  % n
                m0(1) = mean(y./(x.^m0(2)));                  % a
            end
        case 'powerc'    % a*x^n+c;
            % if both x and y are of constant sign:
            if sum(diff(sign(x)))==0 || sum(diff(sign(y)))==0
                m0(3) = log(y(end)/y(1)) / log(x(end)/x(1));  % n
                m0(1) = mean(y./(x.^m0(3)));                  % a
            end
            m0(2) = mean(y);  % c
        case 'powershift'    % a*(x+b)^n;
            % if both x and y are of constant sign:
            m0(2) = mean(x);                              % b
            m0(3) = log(y(end)/y(1)) / log((x(end)+m0(2))/(x(1)+m0(2)));  % n
            m0(1) = mean(y./((x+m0(2)).^m0(3)));                  % a           
        case 'exp'    % a*exp(b*x)
            m0(2) = (log(y(end)/y(1))) / (x(end)-x(1));  % b
            m0(1) = mean(y./exp(m0(2)*x));               % a
        case 'expdiv'    % a*exp(x/b)
            m0(2) = (x(end)-x(1)) / (log(y(end)/y(1)));  % b
            m0(1) = mean(y./exp(x/m0(2)));               % a
        case 'explim'    % a*(1-exp(-x/b))
            m0(1) = max(y);                              % a
            m0(2) = mean(-x./(log(1-y/m0(1))));          % b
        case 'expc'    % a*exp(b*x)+c
            m0(3) = mean(y);                             % c
            m0(2) = (log((y(end)-m0(3))/(y(1)-m0(3)))) / (x(end)-x(1));  % b
            m0(1) = mean((y-m0(3))./exp(m0(2)*x));               % a
        case 'log'    % a*log(b*x)
            m0(1) = (y(end)-y(1))/(log(x(end)/x(1)));    % a
            m0(2) = mean(exp(y/m0(1))./x);               % b            
        case 'logc'    % a*log(x)+b
            m0(1) = (y(end)-y(1))/(log(x(end)/x(1)));    % a
            m0(2) = mean(y - m0(1)*log(x));              % b     
        case {'sin','cos'}    % a*sin(b*x), a*cos(b*x)
            m0(1) = std(y,1)*sqrt(2);  % a
            m0(2) = 50/(x(end)-x(1));  % b
        case {'sinc','cosc'}    % a*sin(b*x)+c, a*cos(b*x)+c
            m0(1) = std(y,1)*sqrt(2);  % a
            m0(2) = 50/(x(end)-x(1));  % b
            m0(3) = mean(y);           % c
        case 'sinphi'    % a*sin(b*x+phi)
            m0(1) = std(y,1)*sqrt(2);  % a
            m0(2) = 50/(x(end)-x(1));  % b
            m0(3) = 1;                 % phi
        case 'sinphic'    % a*sin(b*x+phi)+c
            m0(1) = std(y,1)*sqrt(2);  % a
            m0(2) = 50/(x(end)-x(1));  % b
            m0(3) = mean(y);           % c
            m0(4) = 1;                 % phi
        case 'cngauss'
            m0(1) = (mean(x.^2.*y)/mean(y))^(1/2);   % sigma
        case 'cfgauss'   % a*exp(-(x^2)/(2*sigma^2))
            m0(1) = max(y);      % a
            m0(2) = (mean(x.^2.*y)/mean(y))^(1/2);   % sigma
        case 'ngauss'    % exp(-((x-x_0)^2)/(2*sigma^2))/(2*pi*sigma^2)^(1/2)           
            m0(2) = mean(x.*y)/mean(y);             % x_0
            m0(1) = (mean((x-m0(2)).^2.*y)/mean(y))^(1/2);  % sigma
        case {'fgauss','gauss'}    % a*exp(-((x-x_0)^2)/(2*sigma^2))
            m0(1) = max(y);     % a
            m0(3) = mean(x.*y)/mean(y);             % x_0                   
            m0(2) = (mean((x-m0(3)).^2.*y)/mean(y))^(1/2);  % sigma
    end
    % if an initial guess is 0 or infinity or has imaginary part, set it to 1
    for np=1:length(m0)
        if m0(np)==0 || isinf(m0(np)) || imag(m0(np))~=0
            m0(np) = 1;
        end
    end
end
    

% search for initial guess defined into remfun
remfun = strrep([remfun ';'], ';;', ';'); % adds a final ';'
while strfind(remfun,'='),
    peq=strfind(remfun,'='); peq=peq(1);
    pc=strfind(remfun,';'); pc=pc(1);
    if pc>peq,   % if ';' after '='
        par=remfun(1:(peq-1));
        for i=1:maxm,
            if strcmp(param{i},par),
                m0(i)=eval(remfun((peq+1):(pc-1)));
            end
        end
    else
        error('Ezyfit:ezfit:syntaxError','Invalid syntax');
    end
    remfun=remfun((pc+1):end); % removes the processed i.g. and loops back
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fitting itself
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch f.fitmode,
    case 'lin',
        x__ref = x;
        try
            m=fminsearch(@fitlin, m0);
        catch
            error('Ezyfit:ezfit:fminsearchError','Fit: error during the fminsearch procedure');
        end        
        y_fit = eval(eqml);
        ssr = sum(abs(y_fit-mean(y)).^2);
        sse = sum(abs(y_fit-y).^2);
        f.r=sqrt(ssr/(sse+ssr));
    case 'log',
        x__ref = x;
        try
          m=fminsearch(@fitlog, m0);   % bug fixed here! (v2.51)
        catch
            error('Ezyfit:ezfit:fminsearchError','Fit: error during the fminsearch procedure');
        end
        y_fit = eval(eqml);
        ssr = sum(abs(log(y_fit)-mean(log(y))).^2);
        sse = sum(abs(log(y_fit)-log(y)).^2);
        f.r=sqrt(ssr/(sse+ssr));
    otherwise
        error('Ezyfit:ezfit:unknownMode''Unknown fit mode');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% outputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% fills the output structure:
f.param = param;
f.m = m;
f.m0 = m0;
f.x = x;
f.y = y;
if sum(dy-1) % store the error bars only if defined
    f.dy = dy;
end
if exist('h','var'), f.hdata=h; end  % handle to the data

% stores the fit in the variable 'lastfit' in the 'base' workspace:
assignin('base','lastfit',f);

if strcmp(fp.automakevarfit,'on')
    makevarfit;
end

% ending displays (if no output argument):
if ~nargout
    if strcmp(fp.dispeqmode,'on') % new v2.30
        dispeqfit(f,fp);
    end
    clear f;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of the main function EZFIT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% Nested functions that evaluate the fit for prescribed parameters m(i),
% and return the chi2 (sum of the squared difference between the input
% curve and the fit), in lin or log mode:

    function chi2 = fitlin(m)
        y_fit = eval(eqml);
        chi2 = sum(((y_fit - y).^2)./(dy.^2));
    end

    function chi2 = fitlog(m)
        y_fit = eval(eqml);
        chi2 = sum((log(y_fit)-log(y)).^2);
    end
end