function [f_res,Qunload,Qload,Qext,beta]=qfactor(freq,xx,yy)
% Computes the Q-factor from S11 data and plot the Smith Chart to check the
% results. S11 data needs to be in the detuned short position and they
% should be calibrated
% xx=real(S11)
% yy=imag(S11)
% Choose with the mouse a point at the left of all the intersections
% The method works well for detuned short position, i.e. if the circle
% touches the bigger Smith circle. If it is not the case, only the Q0
% estimate seems to be valid.
% function [f_res,Qunload,Qload,Qext,beta]=qfactor(freq,xx,yy)

% Plot the Smith Chart and the main lines
%Q0 circle
theta=pi/4:0.01:3*pi/4;
Q0circle_xx_up=sqrt(2)*cos(theta);
Q0circle_yy_up=-1+sqrt(2)*sin(theta);

Q0circle_xx_down=sqrt(2)*cos(theta);
Q0circle_yy_down=1-sqrt(2)*sin(theta);

% Qload lines
tt=-1:0.01:0;
yy_up=tt+1;
yy_down=-1-tt;

% Qext circles
alpha=-pi/2:0.01:0;
circle_xx_up=-1+cos(alpha);
circle_yy_up=1+sin(alpha);

circle_xx_down=-1+cos(alpha);
circle_yy_down=-1-sin(alpha);

smithchart(xx+j.*yy);
hold on;
plot(Q0circle_xx_up,Q0circle_yy_up,'r'); 
plot(Q0circle_xx_down,Q0circle_yy_down,'r');
plot(circle_xx_up,circle_yy_up,'g'); 
plot(circle_xx_down,circle_yy_down,'g');
plot(tt,yy_up,'m'); plot(tt,yy_down,'m');
%hold off;

% Finding the intersection point
xmin=-0.7;
[xmin,tmp]=ginput(1);
eps=10^(-2);

abs_diff=abs(xx-xmin);
ii=1;
while((abs_diff(ii)>eps)&&(ii<length(xx)))
    ii=ii+1;
end
jj=ii+200;
while((abs_diff(jj)>eps)&&(jj<length(xx)))
    jj=jj+1;
end
plot(xx(ii),yy(ii),'b+','LineWidth',1)
plot(xx(jj),yy(jj),'b+','LineWidth',1)
aux_xx=xx(ii:jj);
aux_yy=yy(ii:jj);

% Intersection with Q0 circle
diff_circle_up=aux_xx.^2+(aux_yy+1).^2-2;
diff_circle_down=aux_xx.^2+(aux_yy-1).^2-2;
[value_up_Q0,index_up_Q0]=min(abs(diff_circle_up));
[value_down_Q0,index_down_Q0]=min(abs(diff_circle_down));
freq_Q0_up=freq(ii+index_up_Q0-1);
freq_Q0_down=freq(ii+index_down_Q0-1);
xx_Q0_up=xx(ii+index_up_Q0-1);
yy_Q0_up=yy(ii+index_up_Q0-1);
xx_Q0_down=xx(ii+index_down_Q0-1);
yy_Q0_down=yy(ii+index_down_Q0-1);

% Intersection with the Qload lines
diff_lines_up= aux_yy-aux_xx-1;
diff_lines_down= aux_yy+aux_xx+1;
[value_up,index_up_QL]=min(abs(diff_lines_up));
[value_down,index_down_QL]=min(abs(diff_lines_down));
freq_QL_up=freq(ii+index_up_QL-1);
freq_QL_down=freq(ii+index_down_QL-1);
xx_QL_up=xx(ii+index_up_QL-1);
yy_QL_up=yy(ii+index_up_QL-1);
xx_QL_down=xx(ii+index_down_QL-1);
yy_QL_down=yy(ii+index_down_QL-1);

% Intersection with the Qext circles
diff_circle_up2=(aux_xx+1).^2+(aux_yy-1).^2-1;
diff_circle_down2=(aux_xx+1).^2+(aux_yy+1).^2-1;
[value_up_circles,index_up_circles]=min(abs(diff_circle_up2));
[value_down_circles,index_down_circles]=min(abs(diff_circle_down2));
freq_Qext_up=freq(ii+index_up_circles-1);
freq_Qext_down=freq(ii+index_down_circles-1);
xx_circles_up=xx(ii+index_up_circles-1);
yy_circles_up=yy(ii+index_up_circles-1);
xx_circles_down=xx(ii+index_down_circles-1);
yy_circles_down=yy(ii+index_down_circles-1);


% Plot the points for checking
plot(xx_Q0_up,yy_Q0_up,'r*','LineWidth',1)
plot(xx_Q0_down,yy_Q0_down,'r*','LineWidth',1)
plot(xx_QL_up,yy_QL_up,'m*','LineWidth',1)
plot(xx_QL_down,yy_QL_down,'m*','LineWidth',1)
plot(xx_circles_up,yy_circles_up,'g*','LineWidth',1)
plot(xx_circles_down,yy_circles_down,'g*','LineWidth',1)
%hold off

% Finding the final value. Notation according Henke paper.
[min_S11,ind_min]=min(sqrt(xx.^2+yy.^2));
f_res=freq(ind_min);
delta1=(freq_Q0_up-f_res)/freq_Q0_up;
delta2=(freq_Q0_down-f_res)/freq_Q0_down;
delta3=(freq_Qext_up-f_res)/freq_Qext_up;
delta4=(freq_Qext_down-f_res)/freq_Qext_down;
delta5=(freq_QL_up-f_res)/freq_QL_up;
delta6=(freq_QL_down-f_res)/freq_QL_down;

Qunload=-1/(delta1-delta2);
Qload=-1/(delta5-delta6);
Qext=-1/(delta3-delta4);
beta1=Qunload/Qext;
beta2=Qunload/Qload-1;
beta=(beta1+beta2)/2;
plot((beta-1)/(1+beta),0,'rd','LineWidth',0.5)
hold off
