function [H, sigma_H, Cav_LPF, dt, sigma_dt] = Fit_Opgain(fit_ff, fit_abs, fit_phase, Cav, cavname, init_H, model_AS , varargin)
%mn_estimate_Apzt
%  the script to estimate the actuator efficiency of NPRO PZT
[xData, yData] = prepareCurveData( fit_ff*2*pi, fit_abs );

% create string of fitting equation
[z,p,k] = zpkdata(model_AS);
[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);

init_H_TF = tf(init_H,[1/Cav(['Cpole_' cavname])/2/pi 1]);
[z,p,k] = zpkdata(init_H_TF);
[real_zero_index,comp_zero_index] = class_real_comp(z{1});
[real_pole_index,comp_pole_index] = class_real_comp(p{1});
Hrefcav_str = ['(1+dH)*' 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 ')'];
% 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 = 'LAR';

% Fit model to data.
[fitresult, ~] = fit( xData, yData, ft, opts );
H = (1+fitresult.dH)*init_H;
Cav_LPF = (1+fitresult.dH) * init_H_TF;
fit_OLG = Cav_LPF * model_AS;

sigma = confint(fitresult,0.6827);
namelist = coeffnames(fitresult);
sigma_H = abs(init_H*(fitresult.dH-sigma(mnismember(namelist,'dH'))));


[z,p,k] = zpkdata(fit_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);
str_fit = ['unwrap(angle(' str_olg '*exp(-1i*w*dt)))'];
fitting_data_c = fit_abs.*exp(1i*fit_phase*pi/180);
fitting_angle = unwrap(angle(fitting_data_c));

[xData, yData] = prepareCurveData( fit_ff*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 = 'Bisquare';
% Fit model to data.
[fitresult, ~] = fit( xData, yData, ft, opts );
dt = fitresult.dt;
sigma = confint(fitresult,0.6827);
sigma_dt = sigma(1);
if dt <0
    dt = 0;
end
fit_OLG = fit_OLG*tf(1,1,'InputDelay',dt);


if mnismember(varargin,'Plot')
model_c = mnbode_list(fit_OLG,fit_ff,'c');
mnbode(fit_ff,fit_abs,fit_phase,fit_ff,abs(model_c),angle(model_c)*180/pi,...
    'linestyle',{'*','-','-'});
end

end