function [a_pzt, sigma_pzt] = Fit_Apzt( Fitted_Servo_0dB, Act, Act_Sigma, varargin )

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


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

%import data
%Gpzt/(Gaom+1)
for ii = 1:3
data = importdata(['T0831_' num2str(ii) 'G.CSV']);
Freq(['PZT/AOM' num2str(ii)]) = data.data(:,1);
Abs(['PZT/AOM' num2str(ii)]) = 10.^(data.data(:,2)/20);
data = importdata(['T0831_' num2str(ii) 'P.CSV']);
Phase(['PZT/AOM' num2str(ii)]) = data.data(:,2);
end
%Gaom/(Gpzt+1)
%O0904_4: flip the sign of PZT servo.
for ii = 1:4
mnimportdata(['O0904_' num2str(ii) '.txt'],Freq,Abs,Phase,['AOM/PZT' num2str(ii)],'tf');
end

% estimation
A_pzt_list = zeros(1,7);
dA_pzt_list = zeros(1,7);
init_Apzt = -1.3e6;
init_Hrc = 2.8e-6;

% make servo model with CG,FG = 0dB
figure()
for ii = 1:3
    %remove data points with bad coherence
    maxfreq = 2e4;
    [fitting_freq,fitting_abs] = mntrim(0,maxfreq,Freq(['PZT/AOM' num2str(ii)]),Abs(['PZT/AOM' num2str(ii)]));
    [~,fitting_phase] = mntrim(0,maxfreq,Freq(['PZT/AOM' num2str(ii)]),Phase(['PZT/AOM' num2str(ii)]));
    
    % servo TF
    CG = 0;
    FG = -10;
    fast_pol = 1;
    SR560 = -1;
    Servo('TTFSS_PZT') = Fitted_Servo_0dB('TTFSS_PZT')*10^(CG/20)*10^(FG/20)*fast_pol;
    % G_PZT/(G_AOM+1) = C_pztaom fitting
    if ii == 3
        tempdoplot = doplot;
    else
        tempdoplot = '';
    end
    [A_pzt_list(ii),dA_pzt_list(ii)] = estimate_Apzt(fitting_freq,fitting_abs,fitting_phase,init_Hrc,init_Apzt,Servo('TTFSS_PZT'),SR560,Act('Aaom'),tempdoplot);
end



figure()
for ii = 1:3
    %remove data points with bad coherence
    maxfreq = 2e4;
    [fitting_freq,fitting_abs] = mntrim(0,maxfreq,Freq(['AOM/PZT' num2str(ii)]),Abs(['AOM/PZT' num2str(ii)]));
    [~,fitting_phase] = mntrim(0,maxfreq,Freq(['AOM/PZT' num2str(ii)]),Phase(['AOM/PZT' num2str(ii)]));
    % servo TF
    CG = -3;
    FG = -10;
    fast_pol = -1;
    SR560 = -1;
    Servo('TTFSS_PZT') = Fitted_Servo_0dB('TTFSS_PZT')*10^(CG/20)*10^(FG/20)*fast_pol;
    % G_AOM/(G_PZT+1) = C_aompzt fitting
    if ii == 3
        tempdoplot = doplot;
    else
        tempdoplot = '';
    end
    [A_pzt_list(ii+3),dA_pzt_list(ii+3)] = estimate_Apzt2(fitting_freq,fitting_abs,fitting_phase,init_Hrc,init_Apzt,Servo('TTFSS_PZT'),SR560,Act('Aaom'),tempdoplot);
end




figure()
%remove data points with bad coherence
maxfreq = 2e4;
[fitting_freq,fitting_abs] = mntrim(0,maxfreq,Freq('AOM/PZT4'),Abs('AOM/PZT4'));
[~,fitting_phase] = mntrim(0,maxfreq,Freq('AOM/PZT4'),Phase('AOM/PZT4'));
% servo TF
CG = -3;
FG = -10;
fast_pol = 1;
SR560 = -1;
Servo('TTFSS_PZT') = Fitted_Servo_0dB('TTFSS_PZT')*10^(CG/20)*10^(FG/20)*fast_pol;
% G_PZT/(G_AOM+1) = C_pztaom fitting
[A_pzt_list(7),dA_pzt_list(7)] = estimate_Apzt2(fitting_freq,fitting_abs,fitting_phase,init_Hrc,init_Apzt,Servo('TTFSS_PZT'),SR560,Act('Aaom'),doplot);

%result
a_pzt = mean(A_pzt_list);
sigma_pzt = std(A_pzt_list);

end


function [fit_Apzt, sigma_Apzt] = estimate_Apzt(freq, gain, phase, predict_H, predict_Apzt, servo_model,SR560,AOM,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
Hrefcav_str = ['(1+dH)*' num2str(predict_H)];
PZT_str = ['(1+dApzt)*' num2str(predict_Apzt) '*' num2str(SR560)];
AOM_str = num2str(AOM);
[z,p,k] = zpkdata(servo_model);
[real_zero_index,comp_zero_index] = class_real_comp(z{1});
[real_pole_index,comp_pole_index] = class_real_comp(p{1});
servo_str = 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_str '/'...
    '(1+' Hrefcav_str '*' servo_str '*' AOM_str '*exp(-1i*dt*w)))'];

% Set up fittype and options.
ft = fittype( fit_str, 'independent', 'w', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0. 0. 0];
opts.Robust = 'Bisquare';
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
fit_Apzt = (1+fitresult.dApzt)*predict_Apzt;
fit_H = (1+fitresult.dH)*predict_H;
fit_model = fit_H*servo_model*fit_Apzt*SR560 / (1+fit_H*servo_model*AOM*tf(1,1,'InputDelay',fitresult.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','G_{pzt}/(1+G_{aom})','xlim',[min(freq) max(freq)],'legend',{'measured','fitted'});
end

sigma = confint(fitresult,0.6827);
namelist = coeffnames(fitresult);
sigma_Apzt = abs(fit_Apzt-predict_Apzt*(1+sigma(mnismember(namelist,'dApzt'))));

end

function [fit_Apzt, sigma_Apzt] = estimate_Apzt2(freq, gain, phase, predict_H, predict_Apzt, servo_model,SR560,AOM,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
Hrefcav_str = ['(1+dH)*' num2str(predict_H)];
PZT_str = ['(1+dApzt)*' num2str(predict_Apzt)];
AOM_str = num2str(AOM*SR560);
[z,p,k] = zpkdata(servo_model);
[real_zero_index,comp_zero_index] = class_real_comp(z{1});
[real_pole_index,comp_pole_index] = class_real_comp(p{1});
servo_str = 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 '*' AOM_str '/'...
    '(1+' Hrefcav_str '*' servo_str '*' PZT_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 = 'Bisquare';
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
fit_Apzt = (1+fitresult.dApzt)*predict_Apzt;
fit_H = (1+fitresult.dH)*predict_H;
fit_model = fit_H*servo_model*AOM*SR560 / (1+fit_H*servo_model*fit_Apzt);

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

sigma = confint(fitresult,0.6827);
namelist = coeffnames(fitresult);
sigma_Apzt = abs(fit_Apzt-predict_Apzt*(1+sigma(mnismember(namelist,'dApzt'))));

end