function [sys,fitresult,varargout] ...
    = mnTF_fit(ff_fitting_data, abs_fitting_data, ang_fitting_data,...
    model_sys, fitting_zero_list, fitting_pole_list,...
    Sigma,TFname,FitParam,varargin)
%CREATEFIT(FF_FITTING_DATA,ABS_FITTING_DATA)
%  Create a fit.
%
%  Data for 'untitled fit 1' fit:
%      X Input : ff_fitting_dat
%      Y Output: abs_fitting_data
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  See also FIT, CFIT, SFIT.

%  Auto-generated by MATLAB on 15-Jul-2017 18:11:41

% convert model transfer function to string
[z,p,k] = zpkdata(model_sys);
z = z{1}';
p = p{1}';

for ii = 1:length(fitting_zero_list)
    if not(sum(ismember(z(fitting_zero_list),conj(z(fitting_zero_list(ii))))))
        error('complex zero is not chosen as fitting parameter with conjugate pair')
    end
end
for ii = 1:length(fitting_pole_list)
    if not(sum(ismember(p(fitting_pole_list),conj(p(fitting_pole_list(ii))))))
        error('complex pole is not chosen as fitting parameter with conjugate pair')
    end
end



[real_z_index,comp_z_index] = class_real_comp(z);%loacal functions
[real_p_index,comp_p_index] = class_real_comp(p);


tf_str = make_tf_str(z,p,k,fitting_zero_list,fitting_pole_list,real_z_index,comp_z_index,real_p_index,comp_p_index);

% Fit
if mnismember(varargin, 'weight')
    index = mnismember(varargin, 'weight');
    ww = varargin{index+1};
    [xData, yData, weights] = prepareCurveData( ff_fitting_data*2*pi, abs_fitting_data, ww);
else
    [xData, yData] = prepareCurveData(ff_fitting_data*2*pi, abs_fitting_data);
end

% Set up fittype and options.
ft = fittype( tf_str, 'independent', 'w', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Robust = 'Bisquare';
opts.StartPoint = zeros(1,1+length(fitting_zero_list)+length(fitting_pole_list));
opts.Lower = -1*ones(1,1+length(fitting_zero_list)+length(fitting_pole_list));
if mnismember(varargin, 'weight')
    opts.Weights = weights;
end


% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
sigma = confint(fitresult,0.6827);
if sum(mnismember(varargin,'PrintFitResult'))
    fitresult
    gof
end
% Plot fit with data.


fit_coeff_name_list = coeffnames(fitresult);
coeff_zero_list = zeros(1,length(z));
coeff_pole_list = zeros(1,length(p));

for ii = 1:length(fit_coeff_name_list)
    str = fit_coeff_name_list{ii};
    if sum(strfind(str,'dG'))
        fit_G = 1  + getfield(fitresult,str);
        FitParam([TFname '_G']) = k * (fit_G);
        Sigma([TFname '_G']) = abs((getfield(fitresult,str) - sigma(1,ii)) * k);
    elseif sum(strfind(str,'dzero'))
        if sum(strfind(str,'_real'))
            index = str2double(str(11:end));
            coeff_zero_list(index) = coeff_zero_list(index) + (1+getfield(fitresult,str))*real(z(index));
            coeff_zero_list(index+1) = coeff_zero_list(index+1) + (1+getfield(fitresult,str))*real(z(index));
            FitParam([TFname '_zero_real' num2str(index)]) = (1+getfield(fitresult,str))*real(z(index));
            Sigma([TFname '_zero_real' num2str(index)]) = abs((getfield(fitresult,str)-sigma(1,ii))*real(z(index)));
        elseif sum(strfind(str,'_imag'))
            index = str2double(str(11:end));
            coeff_zero_list(index) = coeff_zero_list(index) + 1i * (1+getfield(fitresult,str))*imag(z(index));
            coeff_zero_list(index+1) = coeff_zero_list(index+1) - 1i * (1+getfield(fitresult,str))*imag(z(index));
            FitParam([TFname '_zero_imag' num2str(index)]) = 1i * (1+getfield(fitresult,str))*imag(z(index));
            Sigma([TFname '_zero_imag' num2str(index)]) = 1i*abs((getfield(fitresult,str)-sigma(1,ii))*imag(z(index)));
        else
            index = str2double(str(6:end));
            coeff_zero_list(index) = coeff_zero_list(index) + (1+getfield(fitresult,str))*real(z(index));
            FitParam([TFname '_zero' num2str(index)]) = (1+getfield(fitresult,str))*real(z(index));
            Sigma([TFname '_zero' num2str(index)]) = abs((getfield(fitresult,str)-sigma(1,ii))*real(z(index)));
            
        end
    elseif sum(strfind(str,'dpole'))
        if sum(strfind(str,'_real'))
            index = str2double(str(11:end));
            coeff_pole_list(index) = coeff_pole_list(index) + (1+getfield(fitresult,str))*real(p(index));
            coeff_pole_list(index+1) = coeff_pole_list(index+1) + (1+getfield(fitresult,str))*real(p(index));
            FitParam([TFname '_pole_real' num2str(index)]) = (1+getfield(fitresult,str))*real(p(index));
            Sigma([TFname '_pole_real' num2str(index)]) = abs((getfield(fitresult,str)-sigma(1,ii))*real(p(index)));
        elseif sum(strfind(str,'_imag'))
            index = str2double(str(11:end));
            coeff_pole_list(index) = coeff_pole_list(index) + 1i * (1+getfield(fitresult,str))*imag(p(index));
            coeff_pole_list(index+1) = coeff_pole_list(index+1) - 1i * (1+getfield(fitresult,str))*imag(p(index));
            FitParam([TFname '_pole_imag' num2str(index)]) = 1i * (1+getfield(fitresult,str))*imag(p(index));
            Sigma([TFname '_pole_imag' num2str(index)]) = 1i*abs((getfield(fitresult,str)-sigma(1,ii))*imag(p(index)));
        else
            index = str2double(str(6:end));
            coeff_pole_list(index) = coeff_pole_list(index) + (1+getfield(fitresult,str))*real(p(index));
            FitParam([TFname '_pole' num2str(index)]) = (1+getfield(fitresult,str))*real(p(index));
            Sigma([TFname '_pole' num2str(index)]) = abs((getfield(fitresult,str)-sigma(1,ii))*real(p(index)));
        end
    end
end
sys = zpk([p(fitting_pole_list) coeff_zero_list(find(coeff_zero_list))],[z(fitting_zero_list) coeff_pole_list(find(coeff_pole_list))],fit_G);
plotsys = sys;
if sum(mnismember(varargin,'TimeDelay'))
    fitted_model_sys = minreal(model_sys*sys);
    [f_z,f_p,f_k] = zpkdata(fitted_model_sys);
    fitting_data_c = abs_fitting_data.*exp(1i*ang_fitting_data*pi/180);
    fitting_angle = unwrap(angle(fitting_data_c));
    tf_str = make_tf_str(f_z{1},f_p{1},f_k,[],[],real_z_index,comp_z_index,real_p_index,comp_p_index);
    tf_str = [tf_str(12:end-1) '*exp(-1i*w*dt)'];
    angle_str = ['unwrap(angle(' tf_str '))'];
    
    [xData, yData] = prepareCurveData( ff_fitting_data*2*pi, fitting_angle );

    % Set up fittype and options.
    ft = fittype( angle_str, 'independent', 'w', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.Robust = 'LAR';
    opts.StartPoint = 0;

    % Fit model to data.
    [fitresult_dt, gof_dt] = fit( xData, yData, ft, opts );
    if fitresult_dt.dt<0
        fitresult_dt.dt = 0;
    end
    varargout{1} = fitresult_dt.dt;
    varargout{2} = fitresult_dt;
    plotsys = plotsys * tf(1,1,'InputDelay',varargout{1});
    if sum(mnismember(varargin,'PrintFitResult'))
        fitresult_dt
        gof_dt
    end
end

if sum(mnismember(varargin,'Plot'))
    fitted_model_sys = model_sys*plotsys;
    fitted_model_c = mnbode_list(fitted_model_sys,ff_fitting_data,'complex');
    model_c = mnbode_list(model_sys,ff_fitting_data,'complex');
    figure()
    mnbode(ff_fitting_data,abs_fitting_data,ang_fitting_data,...
        ff_fitting_data,abs(model_c),angle(model_c)*180/pi,...
        ff_fitting_data,abs(fitted_model_c),angle(fitted_model_c)*180/pi,...
        'legend',{'measured','model','fitted'},'linestyle',{'*','-','-'})
    figure()
    semilogx(ff_fitting_data,abs(abs_fitting_data-abs(model_c))./abs_fitting_data,ff_fitting_data,abs(abs_fitting_data-abs(fitted_model_c))./abs_fitting_data)
end
end

function [real_z_index,comp_z_index] = class_real_comp(z)
% the function to classify the real zero/poles and complex zero/poles
% the complex zero/poles list has only one pole/zero of pole/zero conjugate
% pares.
comp_z_index = [];
real_z_index = [];
comp_z = [];
for ii = 1:length(z)
    if isreal(z(ii))
        real_z_index = [real_z_index ii];
    else
        if not(sum(mnismember(comp_z,conj(z(ii)))))
            comp_z = [comp_z z(ii)];
            comp_z_index = [comp_z_index ii];
        end
    end
end
end

function tf_str = make_tf_str(z,p,k,fitting_zero_list,fitting_pole_list,real_z_index,comp_z_index,real_p_index,comp_p_index)

tf_str = ['abs((1+dG)*(' num2str(k) ')*'];
for ii = 1:length(z)
    if sum(mnismember(fitting_zero_list,ii))
        if sum(mnismember(real_z_index,ii))
            tf_str = [tf_str '(1i*w-(' num2str(z(ii)) ')*(1 + dzero' num2str(ii) '))*'];
        elseif sum(mnismember(comp_z_index,ii))
            tf_str = [tf_str '(1i*w-(' num2str(real(z(ii))) ')*(1 + dzero_real' num2str(ii)...
                ')+1i*(' num2str(imag(z(ii))) ')*(1 + dzero_imag' num2str(ii) '))*(1i*w+(' num2str(real(z(ii))) ')*(1 + dzero_real' num2str(ii)...
                ')-1i*(' num2str(imag(z(ii))) ')*(1 + dzero_imag' num2str(ii) '))*'];
        end
    else
        tf_str = [tf_str '(1i*w-(' num2str(z(ii)) '))*'];
    end
end
tf_str = [tf_str(1:end-1) '/'];
for ii = 1:length(p)
    if sum(mnismember(fitting_pole_list,ii))
        if sum(mnismember(real_p_index,ii))
            tf_str = [tf_str '(1i*w-(' num2str(p(ii)) ')*(1 + dpole' num2str(ii) '))/'];
        elseif sum(mnismember(comp_p_index,ii))
            tf_str = [tf_str '(1i*w-(' num2str(real(p(ii))) ')*(1 + dpole_real' num2str(ii)...
                ')+1i*(' num2str(imag(p(ii))) ')*(1 + dpole_imag' num2str(ii) '))/(1i*w+(' num2str(real(p(ii))) ')*(1 + dpole_real' num2str(ii)...
                ')-1i*(' num2str(imag(p(ii))) ')*(1 + dpole_imag' num2str(ii) '))/'];
        end
    else
        tf_str = [tf_str '(1i*w-(' num2str(p(ii)) '))/'];
    end
end
tf_str = [tf_str(1:end-1) ')'];
end