
function [cf,fr]=c_factor(fr,filfile,olfile,calib)
%%

wt=110.*abs((fr./0.0399)./(1+1i.*(fr./1.5915)).^2./(1+1i.*(fr./0.0399)));

%% Filter

figure(4)
dat=read_svk(filfile,1);

fr0=0.0159;
fr1=1.5915;
fr2=31.4535;
fr3=723.4316;
fr4=1.539e3;

fil_mod=(1+i.*(dat(:,2)/fr1)).*(1+i.*(dat(:,2)/fr2))...
      ./(1+i.*(dat(:,2)/fr3))./(1+i.*(dat(:,2)/fr4))./(1+i.*(dat(:,2)/fr0));

ga=dat(:,3)./abs(fil_mod);

l=and(dat(:,5)>0.9999,dat(:,2)<500);
mga1=mean(ga(l));
filg1=mga1.*abs(fil_mod);

subplot(3,1,1)
hold on
loglog(dat(:,2),abs(filg1),'r')
hold off
title(['DC gain: ',num2str(mga1)])

filmod=mga1.*(1+i.*fr/fr1).*(1+i.*fr/fr2)./(1+i.*fr/fr3)./(1+i.*fr/fr4)./(i.*fr/fr0);

%% Openloop

figure(3)
dat=read_svk(olfile,1);

fr0=0.0159;
fr1=1.5915;
fr2=31.4535;
fr3=723.4316;
fr4=1.539e3;

ol_mod=(1+i.*(dat(:,2)/fr1)).*(1+i.*(dat(:,2)/fr2))...
      ./(1+i.*(dat(:,2)/fr3))./(1+i.*(dat(:,2)/fr4))./(1+i.*(dat(:,2)/fr0))...
      ./(dat(:,2)).^2;

%fr0=4.8228;
%ol_mod=1./(1+i.*(dat(:,2)/fr0));

ga=dat(:,3)./abs(ol_mod);

l=and(dat(:,5)>0.99,and(dat(:,2)>10,dat(:,2)<25));
mga=mean(ga(l));
ol=mga.*abs(ol_mod);


mf=min(find(and(dat(:,3)<1,dat(:,5)>0.99)));
ugf=dat(mf,2);
pm=dat(mf,4);

subplot(3,1,1)
hold on
loglog(dat(:,2),abs(ol),'r')
% loglog(fr,abs(ol_mod_mass),fr,abs(ol_mod_pzt),fr,abs(-ol_mod_mass+ol_mod_pzt))
hold off
title(['DC gain: ',num2str(mga),', UGF: ',num2str(ugf),' Hz, Phase mergin: ',num2str(pm),' deg'])
axis([1 1e5 1e-5 1e3])

ol_mod=mga*(1+i.*fr/fr1).*(1+i.*fr/fr2)...
                  ./(1+i.*fr/fr3)./(1+i.*fr/fr4)./(i.*fr/fr0)./fr.^2;

%%

cf=abs(1./wt./filmod.*(1+ol_mod)./calib)';

