function [dG_fss,dG_imc,dt_fss,dt_imc] = FSS_IMC_OLG_fit(ff_fitting_data, abs_fitting_data, ang_fitting_data, FSS_OLG_sys, IMC_OLG_sys ,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_fss,p_fss,k_fss] = zpkdata(FSS_OLG_sys);
z_fss = z_fss{1}';
p_fss = p_fss{1}';

tf_fss_str = make_tf_str(z_fss,p_fss,k_fss,[],[],[],[],[],[]);

[z_imc,p_imc,k_imc] = zpkdata(IMC_OLG_sys);
z_imc = z_imc{1}';
p_imc = p_imc{1}';

tf_imc_str = make_tf_str(z_imc,p_imc,k_imc,[],[],[],[],[],[]);

fit_str = ['unwrap(angle(' tf_fss_str '*(1+(1+dG_imc)*' tf_imc_str '*exp(-1i*w*dt_imc))*exp(-1i*w*dt_fss)))'];

fitting_data_c = abs_fitting_data.*exp(1i*ang_fitting_data*pi/180);
fitting_angle = unwrap(angle(fitting_data_c));
%% Fit
[xData, yData] = prepareCurveData( ff_fitting_data*2*pi, fitting_angle );

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

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
if sum(ismember(varargin,'PrintFitResult'))
    fitresult
    gof
end

% Fit model to data.
if fitresult.dt_imc<0
    fitresult.dt_imc = 0;
end
dt_imc =  fitresult.dt_imc;
dt_fss = fitresult.dt_fss;
dG_imc = 1+fitresult.dG_imc;


%% dG_fss fitting

fit_str = ['abs((1+dG_fss)*' tf_fss_str '*(1+(' num2str(dG_imc) ')*' tf_imc_str '*exp(-1i*w*' num2str(dt_imc) '))*exp(-1i*w*' num2str(dt_fss) '))'];
[xData, yData] = prepareCurveData( ff_fitting_data*2*pi, abs_fitting_data );

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

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
if sum(ismember(varargin,'PrintFitResult'))
    fitresult
    gof
end

dG_fss = 1+fitresult.dG_fss;


if sum(ismember(varargin,'Plot'))
    fitted_model_sys = (dG_fss)*FSS_OLG_sys * (1+(dG_imc)*IMC_OLG_sys*tf(1,1,'InputDelay',dt_imc))*tf(1,1,'InputDelay',dt_fss);
    model_sys = FSS_OLG_sys * (1+IMC_OLG_sys);
    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 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 = ['(' num2str(k) ')*'];
for ii = 1:length(z)
    if sum(ismember(fitting_zero_list,ii))
        if sum(ismember(real_z_index,ii))
            tf_str = [tf_str '(1i*w-(' num2str(z(ii)) ')*(1 + dzero' num2str(ii) '))*'];
        elseif sum(ismember(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(ismember(fitting_pole_list,ii))
        if sum(ismember(real_p_index,ii))
            tf_str = [tf_str '(1i*w-(' num2str(p(ii)) ')*(1 + dpole' num2str(ii) '))/'];
        elseif sum(ismember(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