function [cpole, sigma_cpole] = Fit_Cavpole_IMC( Fitted_Servo_0dB, Act, Cav, Cav_Sigma, varargin )

if mnismember(varargin,'Plot')
    doplot = 'Plot';
else
    doplot = 'noPlot';
end

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

fit_H = zeros(1,10);
sigma_H = zeros(1,10);
fit_cp = zeros(1,10);
sigma_cp = zeros(1,10);

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);
olg_fss = -Cav('LPF_RefCav_0901')*Cav('Mixer_RefCav')*(Servo('TTFSS_PZT')*Act('NPRO_PZT')+Servo('TTFSS_EOM')*Act('BEOM'));

for ii = 1:10
    %importdata
    if ii == 10
        tempdoplot = doplot;
    else
        tempdoplot = '';
    end
    mnimportdata(['T0901_' num2str(ii+40) '.CSV'],Freq,Abs,Phase,['OPGAIN_IMC' num2str(ii)],'tf');
    [~,fitabs] = mntrim(0,2e5,Freq(['OPGAIN_IMC' num2str(ii)]),Abs(['OPGAIN_IMC' num2str(ii)]));
    [fitff,fitphase] = mntrim(0,2e5,Freq(['OPGAIN_IMC' num2str(ii)]),Phase(['OPGAIN_IMC' num2str(ii)]));
    [fit_H(ii),sigma_H(ii),fit_cp(ii),sigma_cp(ii)] =...
        mn_estimate_Himc(fitff,fitabs,fitphase, 1.8e-6, 6e3, Act('Aaom'),minreal(olg_fss),tempdoplot);
end

cpole = mean(fit_cp);
sigma_cpole = std(fit_cp);

sideband_resfreq1 = 11.2486e6;
sideband_resfreq2 = 16.8728e6;
dsideband_resfreq1 = 0.0002e6;
dsideband_resfreq2 = 0.0002e6;
Cav('FSR_IMC') = sideband_resfreq2-sideband_resfreq1;
Cav_Sigma('FSR_IMC') = sqrt(dsideband_resfreq1^2+dsideband_resfreq2^2);
c = 299792000;
Cav('L_IMC') = c/Cav('FSR_IMC')/2;
Cav_Sigma('L_IMC') = Cav('L_IMC')*Cav_Sigma('FSR_IMC')/Cav('FSR_IMC');
Cav('Finnese_IMC') = c/(4*Cav('L_IMC')*cpole);
Cav_Sigma('Finnese_IMC') = Cav('Finnese_IMC')*sqrt((Cav_Sigma('L_IMC')/Cav('L_IMC'))^2+(sigma_cpole/cpole)^2);
F = Cav('Finnese_IMC');
Cav('r_IMC') = (-pi+sqrt(pi^2+4*F^2))/2/F;
Cav('R_IMC') = Cav('r_IMC')^2;
Cav('T_IMC') = 1 - Cav('R_IMC');
Cav_Sigma('r_IMC') = Cav('r_IMC')*(1/(-pi+sqrt(pi^2+4*F^2))*4*F/sqrt(pi^2+4*F^2)-1/F)*Cav_Sigma('Finnese_IMC');
Cav_Sigma('R_IMC') = 2*Cav('R_IMC')*Cav_Sigma('r_IMC');
Cav_Sigma('T_IMC') = Cav_Sigma('R_IMC');
Cav('Loss_IMC') = Cav('T_IMC')-0.0063;
Cav_Sigma('Loss_IMC') = Cav_Sigma('T_IMC');

end

function [fit_H, sigma_H, fit_cavpole, sigma_cavpole] = mn_estimate_Himc(freq, gain, phase, predict_H, predict_cavpole, Aaom, FSS_olg, varargin)
%mn_estimate_Apzt
%  the script to estimate the actuator efficiency of NPRO PZT
[xData, yData] = prepareCurveData( freq*2*pi, gain );

% create string of fitting equation
AOM_str = num2str(Aaom);
Hstr = ['(1+dH)*' num2str(predict_H) '/(1i*w/2/pi/(' num2str(predict_cavpole) '*(1+dcpole))+1)'];
[z,p,k] = zpkdata(FSS_olg/(1+FSS_olg));
[real_zero_index,comp_zero_index] = class_real_comp(z{1});
[real_pole_index,comp_pole_index] = class_real_comp(p{1});
FSS_olg_str = make_tf_str(z{1},p{1},k,[],[],real_zero_index,comp_zero_index,real_pole_index,comp_pole_index);

fit_str = ['abs(' Hstr '*' AOM_str '*' FSS_olg_str ')'];

% Set up fittype and options.
ft = fittype( fit_str, 'independent', 'w', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0. 0. ];
opts.Robust = 'LAR';
% Fit model to data.
[fitresult, ~] = fit( xData, yData, ft, opts );
fit_cavpole = predict_cavpole*(1+fitresult.dcpole);
fit_H = (1+fitresult.dH)*predict_H;
tf_Himc = tf(fit_H,[1/2/pi/fit_cavpole 1]);
fit_model = tf_Himc*Aaom;

sigma = confint(fitresult,0.6827);
namelist = coeffnames(fitresult);
sigma_H = abs(fit_H*(fitresult.dH-sigma(mnismember(namelist,'dH'))));
sigma_cavpole = abs(predict_cavpole*(fitresult.dcpole-sigma(mnismember(namelist,'dcpole'))));


[z,p,k] = zpkdata(fit_model);
[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);
str_fit = ['unwrap(angle(' str_olg '*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.];
opts.Robust = 'LAR';
% Fit model to data.
[fitresult, ~] = fit( xData, yData, ft, opts );
dt = fitresult.dt;
sigma = confint(fitresult,0.6827);
sigma_dt = sigma(1);


fit_model = fit_model*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','LPF of the IMC','xlim',[min(freq) max(freq)],'legend',{'measured','fitted'});
end


end