%==============================================================================================
% Create a Smith chart that can be used to neatly plot a reflection coefficient.
% Smith chart lines calculated using refrence:
% http://www.sss-mag.com/pdf/smith_chart_basics.pdf
%==============================================================================================
% Original by Jackie Cooke
% Updated by James R. Nagel
% Department of Electrical and Computer Engineering
% University of Utah, Salt Lake City, Utah
% Last updated October 1, 2020

function fig = CreateSmithChart(FigNum)

%==============================================================================================
% Initialize figure. 
%==============================================================================================
if nargin == 0
    fig = figure();
else
    fig = figure(FigNum);
    clf;
end

LineWidth = 0.5;
Color = [0.5, 0.5, 0.5];
FontSize = 14;
FontName = 'Helvetica Neue';
FontColor = [0.2, 0.2, 0.2];
hold on

%==============================================================================================
% Plot constant-resistance circles. Remember that resistance circles are all centered at
%   (x0,y0) = [rL/(rL+1),0] ,
% where rL is the normalized load resistance. Likewise, each circle has a radius of
%   r0: 1/(rL+1)
%==============================================================================================

% Define a sampling of representative load resistances.
rL = [0, 0.2, 0.5, 1, 2, 5, 30];

% Calculate circle centers and their radii.
x0 = rL./(rL+1);
r0 = 1./(rL+1);

% Define angular resolution for the circles.
th = linspace(0,2*pi,101); %0:pi/50:2*pi;
cosTH = cos(th);
sinTH = sin(th);

for k = 1:length(rL)
    
    % Calculate (x,y) coordinates for the circles.
    x = r0(k) * cosTH + x0(k);
    y = r0(k) * sinTH;
    
    % Draw the resistance circle.
    p = plot(x,y,'k','LineWidth',LineWidth,'color',Color);
    set(get(get(p,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
    
    % Add text annotation at the far-left end of the resistance circle.
    xT = -r0(k) + x0(k);
    T = text(xT,0,num2str(rL(k)));
    set(T,'Color',FontColor,'FontSize',FontSize,'FontName',FontName);
    set(T,'horizontalalignment','right');
end

% Add a final text annotation for rL = infinity and place it at the right edge of the unit
% circule.
T = text(1,0,'\infty');
set(T,'Color',FontColor,'FontSize',FontSize,'FontName',FontName);
set(T,'horizontalalignment','left');

%==============================================================================================
% Plot constant-reactance circles. Remember that reactance circles are all centered at
%   (x0,y0) = [1,1/xL] ,
% where xL is the normalized load reactance. Likewise, each circle has a radius of
%   r0: 1/xL
% Note that we are only interested in the part of the reactance circle that falls within the
% unit circle. So really, each reactance circle is only going to be a "reactance arc."
%==============================================================================================

% Begin with inductive reactance circles. These are all defined by positive xL and lie above
% the x-axis.  
xL = [0, 0.2, 0.5, 1, 2, 5, 30];

% Define how many angular samples to include with each reactance arc.
Nth = 101;

% Calculate circle centers and their radii.
x0 = 1;
y0 = 1./xL;
r0 = 1./xL;

% Calculate angular span of the reactance circle that falls within the unit circle. Note that
% the derivation for this basically just intersecting two circles and applying some clever
% trigonometry. 
phi = 2*atan(1./r0);

% Now convert the angular span into start/stop angles for the reactance circle.
thetaStart = 3*pi/2 - phi;
thetaStop  = 3*pi/2;

for k = 1:length(xL)
  
    % Build an array of angular samples to plot the arc.
    TH = linspace(thetaStart(k),thetaStop,Nth);
    
    % Calculate (x,y) coordinates for the circular arc.
    x = r0(k)*cos(TH) + x0;
    y = r0(k)*sin(TH) + y0(k);
    p = plot(x,y,'-','LineWidth',LineWidth,'color',Color);
    set(get(get(p,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
    
    % Plot the corresponding reactance arc. By symmetry, this will be the same circle, but
    % with the y-coordinates negated.
    p = plot(x,-y,'-','LineWidth',LineWidth,'color',Color);
    set(get(get(p,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
    
    % Insert text annotation for the positive reactance circle. Place each annotation at the
    % edge of the unit circle and tweak alignment for neatness.
    xTip = x(1);
    yTip = y(1);
    ArcLabel = ['+j' num2str(xL(k))];
    
    T = text(xTip,yTip,ArcLabel);
    set(T,'Color',FontColor,'FontSize',FontSize,'FontName',FontName);
    set(T,'VerticalAlignment','bottom');
    if xTip < 0
        set(T,'horizontalalignment','right');
    elseif xTip == 0
        set(T,'horizontalalignment','center');
    else
        set(T,'horizontalalignment','left');
    end
   
    % Insert text annotation for the negative reactance circle. Place each annotation at the
    % edge of the unit circle and tweak alignment for neatness.
    xTip = x(1);
    yTip = -y(1);
    ArcLabel = ['-j' num2str(xL(k))];
    
    T = text(xTip,yTip,ArcLabel);
    set(T,'Color',FontColor,'FontSize',FontSize,'FontName',FontName);
    set(T,'VerticalAlignment','top');
    if xTip < 0
        set(T,'horizontalalignment','right');
    elseif xTip == 0
        set(T,'horizontalalignment','center');
    else
        set(T,'horizontalalignment','left');
    end
    
end

% Polish up the axes.
axis([-1 1 -1 1]);
axis equal
axis off
hold off

end