function [fitting_pzt, fitting_eom, fitservo_0dB_pzt, fitservo_0dB_eom, fitservo_0dB_slow] = Fit_TTFSS( TTFSS_Elec, AA_AI, FitParam, Sigma, varargin )

Model = containers.Map;
Freq = containers.Map;
Abs = containers.Map;
Phase = containers.Map;
temp_FitParam = containers.Map;
temp_Sigma = containers.Map;


% plot or not
if mnismember(varargin,'Plot')
    doplot = true;
else
    doplot = false;
end

%import fitting datas
data = importdata('T0712_1G.CSV');
Freq('TTFSS_EOM') = data.data(:,1);
Abs('TTFSS_EOM') = 10.^(data.data(:,2)/20);
data = importdata('T0712_1P.CSV');
Phase('TTFSS_EOM') = data.data(:,2);
data = importdata('T0712_2G.CSV');
Freq('TTFSS_PZT') = data.data(:,1);
Abs('TTFSS_PZT') = 10.^(data.data(:,2)/20);
data = importdata('T0712_2P.CSV');
Phase('TTFSS_PZT') = data.data(:,2);

% make model of servos
Model('TTFSS_PZT') = minreal(TTFSS_Elec('common_gain1') * TTFSS_Elec('common_gain2') ...
    * TTFSS_Elec('fast_gain1') * TTFSS_Elec('fast_LPF_10k') * TTFSS_Elec('fast_LPF_34k')...
    * TTFSS_Elec('fast_notch') * TTFSS_Elec('fast_LPF_out'));
Model('TTFSS_EOM') = minreal(TTFSS_Elec('common_gain1') * TTFSS_Elec('common_gain2') ...
    * (TTFSS_Elec('EOM_BPF1') * TTFSS_Elec('EOM_BPF2') * TTFSS_Elec('EOM_HPF') + TTFSS_Elec('EOM_BPF3')) * TTFSS_Elec('EOM_HV'));

FB_filter = mnzpk([],-[100 0],10^(-93/20));
Model('TTFSS_SLOW') = TTFSS_Elec('common_gain1') * TTFSS_Elec('common_gain2') ...
    * TTFSS_Elec('fast_gain1') * TTFSS_Elec('fast_LPF_10k') * TTFSS_Elec('fast_LPF_34k')...
    * TTFSS_Elec('fast_notch') * (TTFSS_Elec('Interface_PZT_gain_p') - (-TTFSS_Elec('Interface_PZT_gain_n')))...
    * AA_AI('Analog_AA') * 2^16/40 * AA_AI('Digital_AA') * FB_filter * AA_AI('Analog_AI') * 20/2^16 * AA_AI('Digital_AI')...
    * 2 * TTFSS_Elec('Interface_SLOW_LPF') * (TTFSS_Elec('SLOW_output_p') - (-TTFSS_Elec('SLOW_output_n')));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PZT path
% remove bad coherence data
[~,abs_fitting_data] = mntrim(0,5e4,Freq('TTFSS_PZT'),Abs('TTFSS_PZT'));
[ff_fitting_data,ang_fitting_data] = mntrim(0,5e4,Freq('TTFSS_PZT'),Phase('TTFSS_PZT'));

% define gains and sign for fitting
FG = -4;
CG = -4;

% devide the fitting data by valuable gains.
abs_fitting_data = abs_fitting_data/10^(FG/20)/10^(CG/20);

% Fitting
[fitting_pzt,~,~] = ...
    TF_Fit(ff_fitting_data,abs_fitting_data,ang_fitting_data,...
    Model('TTFSS_PZT'),[],[1 2],temp_Sigma,'TTFSS_PZT',temp_FitParam,'TimeDelay');
FitParam('TTFSS_PZT_10kHz_pole') = temp_FitParam('TTFSS_PZT_pole1');
Sigma('TTFSS_PZT_10kHz_pole') = temp_Sigma('TTFSS_PZT_pole1');
FitParam('TTFSS_PZT_34kHz_pole') = temp_FitParam('TTFSS_PZT_pole2');
Sigma('TTFSS_PZT_34kHz_pole') = temp_Sigma('TTFSS_PZT_pole2');
FitParam('TTFSS_PZT_G') = temp_FitParam('TTFSS_PZT_G');
Sigma('TTFSS_PZT_G') = temp_Sigma('TTFSS_PZT_G'); 

fitservo_0dB_pzt = minreal(Model('TTFSS_PZT') * fitting_pzt);

if doplot
    close all
    freq = 10.^(0:0.02:6);
    model_c = mnbode_list(fitservo_0dB_pzt,freq,'c');
    mnbode(ff_fitting_data,abs_fitting_data,angle(-exp(1i*pi/180*ang_fitting_data))*180/pi,...
        freq,abs(model_c),angle(model_c)*180/pi,...
        'ylim',[1e-7 5e2],'title','TTFSS OUT1 to FAST (CG,FG = -4dB, fast pol:+)','linestyle',{'*','-'},'legend',{'measured','fitted'})
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EOM path
% remove bad coherence data
[~,abs_fitting_data] = mntrim(0,3e4,Freq('TTFSS_EOM'),Abs('TTFSS_EOM'));
[ff_fitting_data,ang_fitting_data] = mntrim(0,3e4,Freq('TTFSS_EOM'),Phase('TTFSS_EOM'));

% devide the fitting data by valuable gains.
abs_fitting_data = abs_fitting_data/10^(CG/20);

[fitting_eom,~] = ...
    TF_Fit(ff_fitting_data,abs_fitting_data,ang_fitting_data,...
    minreal(Model('TTFSS_EOM')),[],[1 3 10],temp_Sigma,'TTFSS_EOM',temp_FitParam);
FitParam('TTFSS_EOM_2kHz_pole') = temp_FitParam('TTFSS_EOM_pole1');
Sigma('TTFSS_EOM_2kHz_pole') = temp_Sigma('TTFSS_EOM_pole1');
FitParam('TTFSS_EOM_1kHz_pole') = temp_FitParam('TTFSS_EOM_pole3');
Sigma('TTFSS_EOM_1kHz_pole') = temp_Sigma('TTFSS_EOM_pole3'); 
FitParam('TTFSS_EOM_33kHz_pole') = temp_FitParam('TTFSS_EOM_pole10');
Sigma('TTFSS_EOM_33kHz_pole') = temp_Sigma('TTFSS_EOM_pole10'); 
FitParam('TTFSS_EOM_G') = temp_FitParam('TTFSS_EOM_G');
Sigma('TTFSS_EOM_G') = temp_Sigma('TTFSS_EOM_G'); 

fitservo_0dB_eom = minreal(Model('TTFSS_EOM') * fitting_eom);

if doplot
    close all
    freq = 10.^(0:0.02:7);
    model_c = mnbode_list(fitservo_0dB_eom,freq,'c');
    mnbode(ff_fitting_data,abs_fitting_data,ang_fitting_data,...
        freq,abs(model_c),angle(model_c)*180/pi,...
        'ylim',[1e-4 1e6],'xlim',[1 1e7],...
        'title','TTFSS OUT1 to EOM (CG,FG = -4dB)','linestyle',{'*','-'},'legend',{'measured','fitted'})
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SLOW path
fitservo_0dB_slow = minreal(Model('TTFSS_SLOW'));


end