%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%             Thank you for using ClarkWare Products!!!!               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SMITH    Smith Chart.
%          SMITH(Zl,Zo) draws a Smith chart of the Complex impedances
%          in Zl, normalized to the real valued characteristic impedance Zo.
%          Multiple plots are obtained in different colours etc if Zl 
%          is a matrix, ie its columns  contain sets of complex values. 
%          Smith may be subplotted with fair deal of success. (Text size
%          doesn't change and is a bit messy with subplot (221) -- but 
%          liveable! ) Obeys all usual conventions.
%
 
% Originally based on some Matlab code by 
%      Antony-Dean McKechnie & Neville Wilken in their final year project
%      of their BSc(Eng) at Wits. (Which in turn was based on some Pascal
%      code of mine :-)
% But rather heavily rewritten :-) by 
%      Alan Robert Clark, PrEng.
%      Dept Elec Eng
%      P.O.Wits
%      2050 South Africa
%      clark@odie.wits.ee.ac.za
%      15 December 1992.
% Added: subplot capability.
%        multiple plot capability (re-defaulted).
%        'Polar'ized axes
%        combined multi stage process into one simple function call.
% Removed: non standard calling conventions.
%          Bugs, bugs, bugs :-) :-) Probably should fall under Added too!

function smith (Zl, Zo);
% Set up Constant Resistance Circles; 
% and Constant reactance circles -- but crop at Outer Circle.
r = [1 0.826 0.665 0.500 0.334 0.162];
axis([-1 1 -1 1]);
Re = [1,2,5,5];
Im = [0.2 0.5 1 2];
for i = 1:4
 p = (0:(Re(i)/100):Re(i))';
 q = p-p+Im(i);
  Icirc1 = p+j*q;
  Icirc2 = p-j*q;
  rho1=(Icirc1-1)./(Icirc1+1);
  rho2=(Icirc2-1)./(Icirc2+1);
  rhomag = abs(rho1);
  Circm(:,i)= rhomag;
  Circp = angle(rho1);
  Circp1(:,i) = Circp;
  Circp = angle(rho2);
  Circp2(:,i) = Circp;
end,

inc = 2*pi/100;
theta = 0:inc:pi;
phi = 0:2*inc:2*pi;
RV = phi-phi+0.329;
for i=1:6
 R1 = 2*r(:,i)*sin(theta+pi/2);
 y1 = R1.*sin(theta);
 x1 = R1.*cos(theta);
 y2 = y1;
 x2 = x1+(1-2*r(i));
 R = sqrt(x2.^2+y2.^2);
 R2(:,i) = R';
 theta1 = atan((y2./x2));
for p=1:length(x2) % Crop at 1
 if x2(p) < 0, theta1(p)=theta1(p)-pi; end;
end,
 theta2(:,i) = theta1';
end

% Keep polar from doing anything silly
hold on
polar(theta2,R2,':c1',Circp1,Circm,':c1',....
Circp2,Circm,':c1',phi,RV,':c1');

% Fill in the Real axis..
x=[-pi 0];
y=[1 1];
polar (x,y,'w:');

% Do a bit of fiddling to get the scale positioned..
text(-pi,1.1,'0');
text(-pi,0.66,'.2');
text(-pi,0.33,'.5');
text(0,0,'1');
text(0,0.33,'2');
text(0,0.66,'5');
%text(0,1.1,'ì'); % Works fine on screen but PCGPP falls over and dies!!!
text(0,1.05,'inf'); % so we have to make do with a poor man's infinity!!
text(5.2*pi/6,1.1,'.2');
text(-5.2*pi/6,1.15,'.2');
text(2.2*pi/3,1.1,'.5');
text(-2.2*pi/3,1.15,'.5');
text(pi/2,1.05,'1');
text(-pi/2,1.1,'1');
text(pi/3.2,1.05,'2');
text(-pi/3.2,1.1,'2');

% OK now go ahead and plot the impedances, -- first normalize
rho=(Zl-Zo)./(Zl+Zo);
rhomag = abs(rho);
rhoph = angle(rho);
polar(rhoph,rhomag);
hold off % return to usual state 

% Need to address a serious bug in the reversion to rectangular coords 
% .. if the plot is held, and PLOT is used, the FIRST point is haywire--
% causing GPP to hang the machine etc. SECOND point is fine.
% Hence :
hold on; plot(0,0,'i'); hold off;
% which sorts it out beautifully.....



function [rads] = angle (comp_val)
% Can you believe that ANGLE is not defined????
% Angle of a complex value returns radians.....
% Alan Robert Clark
rads = atan2 (imag(comp_val),real(comp_val));