function [mixer_pole, sigma_mixer_pole, dt_eom, sigma_dt_eom, mixer_tf]...
    = Fit_MixerPole(TTFSS_Elec, Fitted_Servo_0dB, Cav, Act, varargin)
%UNTITLED9 Summary of this function goes here
%   Detailed explanation goes here
% plot or not
if mnismember(varargin,'Plot')
    doplot = 'Plot';
else
    doplot = 'noPlot';
end

Freq = containers.Map;
Abs = containers.Map;
Phase = containers.Map;
Servo = containers.Map;


fast_pol = -1;
CG = 30;
FG = 30;
Servo('TTFSS_PZT') = Fitted_Servo_0dB('TTFSS_PZT')*10^(CG/20)*10^(FG/20)*fast_pol;
Servo('TTFSS_EOM') = Fitted_Servo_0dB('TTFSS_EOM')*10^(CG/20);

model_olg = -Cav('LPF_RefCav_0901')*(Servo('TTFSS_PZT')*Act('NPRO_PZT')+Servo('TTFSS_EOM')*Act('BEOM'));


fit_mp = zeros(1,17);
sigma_mp = zeros(1,17);
fit_dt = zeros(1,17);
sigma_dt = zeros(1,17);
for ii = 1:5
    %import data
    mnimportdata(['O0901_' num2str(ii) '.CSV'],Freq,Abs,Phase,['OLG_FSS' num2str(ii)],'tf');
    [~,fitabs] = mntrim(0,2e6,Freq(['OLG_FSS' num2str(ii)]),Abs(['OLG_FSS' num2str(ii)]));
    [fitff,fitphase] = mntrim(0,2e6,Freq(['OLG_FSS' num2str(ii)]),Phase(['OLG_FSS' num2str(ii)]));
    
    [fit_mp(ii),sigma_mp(ii),fit_dt(ii),sigma_dt(ii)] = ...
        mn_estimate_OLGdt(fitff,fitabs*-TTFSS_Elec('common_gain1'),fitphase,...
        1.6e6,model_olg,doplot);
end
for ii = 4:15
    %import data
    mnimportdata(['O0831_' num2str(ii) '.CSV'],Freq,Abs,Phase,['OLG_FSS' num2str(ii+2)],'tf');
    [~,fitabs] = mntrim(0,2e6,Freq(['OLG_FSS' num2str(ii+2)]),Abs(['OLG_FSS' num2str(ii+2)]));
    [fitff,fitphase] = mntrim(0,2e6,Freq(['OLG_FSS' num2str(ii+2)]),Phase(['OLG_FSS' num2str(ii+2)]));
    
    [fit_mp(ii+2),sigma_mp(ii+2),fit_dt(ii+2),sigma_dt(ii+2)] = ...
        mn_estimate_OLGdt(fitff,fitabs*-TTFSS_Elec('common_gain1'),fitphase,...
        1.6e6,model_olg,doplot);

end
mixer_pole = mean(abs(fit_mp));
sigma_mixer_pole = std(abs(fit_mp));
dt_eom = mean(abs(fit_dt));
sigma_dt_eom = std(abs(fit_dt));
mixer_tf = tf(1,[1/2/pi/mixer_pole 1])^2;
end

function [fit_mixerpole,sigma_mixerpole,dt,sigma_dt] = mn_estimate_OLGdt(freq, gain, phase, pre_mixerpole, model_olg,varargin)
%mn_estimate_Apzt

% create string of fitting equation

mixer_str = ['1/(1i*w/(' num2str(pre_mixerpole*2*pi) '*(1+dmpole))+1)^2'];
[z,p,k] = zpkdata(model_olg);
[real_zero_index,comp_zero_index] = class_real_comp(z{1});
[real_pole_index,comp_pole_index] = class_real_comp(p{1});
str_olg = make_tf_str(z{1},p{1},k,[],[],real_zero_index,comp_zero_index,real_pole_index,comp_pole_index);

fit_str = [str_olg '*' mixer_str];


str_fit = ['unwrap(angle(' fit_str '*exp(-1i*w*dt)))'];
fitting_data_c = gain.*exp(1i*phase*pi/180);
fitting_angle = unwrap(angle(fitting_data_c));

[xData, yData] = prepareCurveData( freq*2*pi, fitting_angle );
% Set up fittype and options.
ft = fittype( str_fit, 'independent', 'w', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0. 0.];
opts.Robust = 'Bisquare';
opts.Lower = [-1 -1];
opts.Upper = [1 1];
% Fit model to data.
[fitresult, ~] = fit( xData, yData, ft, opts );
dt = fitresult.dt;
fit_mixerpole = pre_mixerpole*(1+fitresult.dmpole);
fit_mixer = zpk([],-[fit_mixerpole fit_mixerpole]*2*pi,fit_mixerpole^2*2^2*pi^2);

sigma = confint(fitresult,0.6827);
namelist = coeffnames(fitresult);
sigma_mixerpole = abs(fitresult.dmpole-sigma(mnismember(namelist,'dmpole')))*pre_mixerpole;
sigma_dt = sigma(mnismember(namelist,'dt'));

fit_model = fit_mixer*model_olg*tf(1,1,'InputDelay',dt);
if mnismember(varargin,'Plot')
model_c = mnbode_list(fit_model,freq,'c');
mnbode(freq,gain,phase,freq,abs(model_c),angle(model_c)*180/pi,...
    'linestyle',{'*','-'},'title','OLG fitting');
end

end