function [act_eom, sigma_act_eom, tf_eom] = Fit_Aeom(FitServo_0dB, Act, Act_Sigma, Cav, 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;

CG = 30;
FG = 30;
fast_pol = -1;
init_Aeom = -0.00255;

% importdata
mnimportdata('O0901_6.txt',Freq,Abs,Phase,'CROSS_FSS','tf');
Servo('TTFSS_PZT') = FitServo_0dB('TTFSS_PZT')*10^(CG/20)*10^(FG/20)*fast_pol;
Servo('TTFSS_EOM') = FitServo_0dB('TTFSS_EOM')*10^(CG/20);
[act_eom,sigma] = mn_estimate_Aeom(Freq('CROSS_FSS'),Abs('CROSS_FSS'),Phase('CROSS_FSS'),Cav('LPF_RefCav_0901'),...
    init_Aeom, Servo('TTFSS_PZT'),Servo('TTFSS_EOM'),Act('NPRO_PZT'),doplot);
sigma_act_eom = act_eom * sqrt((Act_Sigma('NPRO_PZT')/Act('NPRO_PZT'))^2+(sigma/act_eom)^2);
tf_eom = act_eom*tf([1 0],[1/1e9 1]);
saveas(gcf,'figure/RefCav_opgain.pdf')
end

function [fit_Aeom, sigma_Aeom] = mn_estimate_Aeom(freq, gain, phase, Hpc, predict_Aeom, servo_model_pzt,servo_model_eom,Apzt,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
[z,p,k] = zpkdata(Hpc);
[real_zero_index,comp_zero_index] = class_real_comp(z{1});
[real_pole_index,comp_pole_index] = class_real_comp(p{1});
Hrefcav_str = make_tf_str(z{1},p{1},k,[],[],real_zero_index,comp_zero_index,real_pole_index,comp_pole_index);

EOM_str = ['(1+dAeom)*1i*w*' num2str(predict_Aeom)];
PZT_str = num2str(Apzt);
[z,p,k] = zpkdata(servo_model_pzt);
[real_zero_index,comp_zero_index] = class_real_comp(z{1});
[real_pole_index,comp_pole_index] = class_real_comp(p{1});
servo_str_pzt = make_tf_str(z{1},p{1},k,[],[],real_zero_index,comp_zero_index,real_pole_index,comp_pole_index);
[z,p,k] = zpkdata(servo_model_eom);
[real_zero_index,comp_zero_index] = class_real_comp(z{1});
[real_pole_index,comp_pole_index] = class_real_comp(p{1});
servo_str_eom = make_tf_str(z{1},p{1},k,[],[],real_zero_index,comp_zero_index,real_pole_index,comp_pole_index);

fit_str = ['abs(' Hrefcav_str '*' servo_str_pzt '*' PZT_str '/'...
    '(1+' Hrefcav_str '*' servo_str_eom '*' EOM_str '))'];

% Set up fittype and options.
ft = fittype( fit_str, 'independent', 'w', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0.];
opts.Robust = 'Bisquare';
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
fit_Aeom = (1+fitresult.dAeom)*predict_Aeom;
tf_Aeom = fit_Aeom*tf([1 0],[1/1e9 1]);
fit_model = Hpc*servo_model_pzt*Apzt / (1+Hpc*servo_model_eom*tf_Aeom);

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','G_{pzt}/(1+G_{eom})','legend',{'measured','fit'},'xlim',[min(freq) max(freq)]);
end

sigma = confint(fitresult,0.6827);
namelist = coeffnames(fitresult);
sigma_Aeom = abs(fitresult.dAeom-sigma(mnismember(namelist,'dAeom')))*predict_Aeom;

end