Plot Smith Chart Data in 3-D with MATLAB

MATLAB is a wonderful tool for the capture, analysis, and display of data.  And their Home version provides an inexpensive way for the casual user (such as myself) to utilize its analytical power.

Recently, while analyzing Antenna Tuner power dissipation versus load impedance, I wanted some way to visually represent this loss versus load impedance. 

Because the Smith Chart conveniently represents an infinite range of impedances within a unit-circle (R from 0 to +infinity, and X from -infinity to +infinity), it would provide the ideal co-ordinate system onto which I could map my loss values.

In doing so, I would essentially be adding a third axis to the Smith Chart -- an axis that represents the magnitude of a value that is a function of the Smith Chart's R and X coordinates.  I.e. V = f(R, X), where V is the value to be plotted and f is the function in terms of Resistance and Reactance.

 

How to display three axis on a 2-D plot?  Use MATLAB, of course.

I first attempted to use MATLAB's smithchart() function (in their RF Toolbox) with their surfc() function, but without success.

Fortunately, Dick Benson (W1QG) had written his own Smith Chart function (having been frustrated with the one provided in the RF Toolbox), and he showed me how his function could be used for my requirements

Below is an example.  It is a plot displaying power loss in low-pass L-networks, given component Q (and calculating the L and C values for a 1:1 match) for SWRs from 1:1 to 10:1.  Thus, the outer circumference of the "disk" of calculated power-loss corresponds to L-network load impedances lying on a Smith Chart's 10 to 1 SWR circle.  (Note that low loss is blue, and high loss is red.)

Contour lines (that match the index "tick" marks on the right-hand color bar scale) are also displayed.

This plot can be rotated in 3-D space:

Dick's enhanced Smith Chart function (smith_rab_v2.m) is available on the MATLAB file-central website in his S-Parameter Utilities zip file.  This file includes a documentation PDF.

But Dick's documentation only describes normal "2-D" plots.  To plot in 3-D, we need to use another MATLAB function, surfc():

surfc

Contour plot under a 3-D shaded surface plot

collapse all in page

Syntax

surfc(Z)
surfc(Z,C)
surfc(X,Y,Z)
surfc(X,Y,Z,C)
surfc(...,'PropertyName',PropertyValue)
surfc(axes_handles,...)
h = surfc(...)

Description

surfc(Z) creates a contour plot under the three-dimensional shaded surface from the z components in matrix Z, using x = 1:n and y = 1:m, where[m,n] = size(Z). The height, Z, is a single-valued function defined over a geometrically rectangular grid. Z specifies the color data, as well as surface height, so color is proportional to surface height.

surfc(Z,C) plots the height of Z, a single-valued function defined over a geometrically rectangular grid, and uses matrix C, assumed to be the same size as Z, to color the surface.

surfc(X,Y,Z) uses Z for the color data and surface height. X and Y are vectors or matrices defining the x and y components of a surface. If X and Yare vectors, length(X) = n and length(Y) = m, where [m,n] = size(Z). In this case, the vertices of the surface faces are (X(j), Y(i), Z(i,j)) triples. To create X and Y matrices for arbitrary domains, use the meshgrid function.

surfc(X,Y,Z,C) uses C to define color. MATLAB® performs a linear transformation on this data to obtain colors from the current colormap.

surfc(...,'PropertyName',PropertyValue) specifies surface properties along with the data.

surfc(axes_handles,...) plots into the axes with handle axes_handle instead of the current axes (gca).

h = surfc(...) returns handles to a chart surface and a contour object.

3-D Plotting using smith_rab_v2.m:

Below is the code to create a Smith Chart plot of data in 3-D.  Note that it invokes surfc in the form surfc(X,Y,Z), per the description, above.

% Plot 3-d smith chart

% First, set Smith Chart parameters for smith_rab_v2.m,
% per its documentation:

LW=2;  % Line Width
SP.Rvalues=[0 10 25 50 100 250];   % SP.*= Smith Parameters
SP.Xvalues=[10 25 50 100 200 500];
SP.Zo=50;
SP.Nseg=60;
SP.LW  = 1;                 % LW= line width for Smith coordinates
SP.swr_circles = [];        % set to [] for no SWR circles
SP.LW_swr = 2;              % LW= line width for SWR Circles

SP.colors.grid   = [1 1 0]*0.5;       % coordinate lines
SP.colors.fill   = [0,0,0];           % background
SP.colors.inner_text = [1 1 0];       % inner labels
SP.colors.outer_text = [1 1 0];       % outer labels
SP.colors.swr    = [1 0 0];           % swr circle color
SP.colors.Q      = [0 0 1];           % Q contour color

SP.Q_contours    = [];
SP.Q_pts         = [];  % number of points in a Q contour
SP.LW_Q          = 1;

my_color_choice = [0.5 0.5 0.5];

% And now, PLOT!
handles.figure(4)     = figure('NumberTitle','Off',...
                        'Name', ' Match Space',...
                        'color',my_color_choice,...
                        'Position',[380 300 756 695]);
handles.Axes_Smith(4) = axes;
handles.smith(4)      = smith_rab_v2(handles.Axes_Smith(4),SP);
set(handles.Axes_Smith(4),'visible','on',...
                          'Xcolor',my_color_choice,...
                          'Ycolor',my_color_choice,...
                          'Zcolor',SP.colors.outer_text);
zlabel(['Gp Power Loss in PERCENT']);
title(handles.Axes_Smith(4),['Gp Percent Power Loss'],'Color',[1 1 0]);
colormap jet;

handles.cb(1)=colorbar;
set(handles.cb(1),'Color', SP.colors.outer_text,...
                  'FontWeight','Bold');
set(handles.cb(1).Label,'String','Percent Power Loss',...
                  'Color',SP.colors.outer_text,...
                  'FontWeight','Bold');

handles.surf(1,:) = surfc(real_gamma,imag_gamma,Percent_Loss);
set(handles.surf(1,1),'Edgecolor','none',...
                      'FaceAlpha',0.5,...
                      'LineStyle','none');
set(handles.surf(1,2),'Linewidth',2);  % nice fat contours lines!

sinfo = {['Ql : ',num2str(Q_L)],...
         ['Qc : ',num2str(Q_C)]};

ht=text(-1,1,15,sinfo);
set(ht,'Color',[0 1 0],'FontSize',12)

Notes on the plotting code, above:

1.  The code begins by customizing the Smith Chart plot parameters for smith_rab_v2.m (refer to the documentation PDF contained in the  S-Parameter Utilities download).

2.  The actual 3-D plot is made using MATLAB's surfc.  In the script below, I call surfc in the form of surfc(X,Y,Z):

handles.surf(1,:) = surfc(real_gamma,imag_gamma,Percent_Loss);

where, given the function's input variable's X, Y, Z (i.e. surfc(X,Y,Z)):  

• X = real_gamma (i.e. an array containing the Cartesian x-coordinates of gamma (in this particular example, the array is a 2-D matrix)).
• Y = imag_gamma (i.e. an array containing the Cartesian y-coordinates of gamma (in this particular example, the array is a 2-D matrix)).
• Z = the percent value of power loss for each x,y pair.  This data is also an array, and in this particular example it is a 2-D matrix whose dimensions are the same as the real_gamma and the imag_gamma matrices.

(For a more detailed description of these three matrices and how they can be created, see the comments in the comment section of this blog post).

3.  Note that to be able to view the Smith Chart axes that lie beneath the 3-D surface overlay, the plotted surface must be transparent.  Transparency is set via the 'FaceAlpha' parameter in the line:

set(handles.surf(1,1),'Edgecolor','none',...
                      'FaceAlpha',0.5,...
                      'LineStyle','none');

4.  The parameter 'colormap' determines what colors that data values should be mapped into.  The map I've used is the predefined 'jet' map, which colors the data from blue for the lowest values of data to red for the highest values.

5.  Hopefully the code is self-explanatory.  But if something is puzzling, note that MATLAB provides excellent documentation of their functions, which you can either google (e.g. google 'Matlab surfc') or invoke by using the 'doc' command in MATLAB's workspace (e.g. type 'doc surfc' in the Workspace command line).


Applying 3-D Plotting:

Below is the complete MATLAB script I used to generate the power-dissipation data I've displayed, above.  This code analyzes both the series-L/parallel-C and parallel-C/series-L (LsCp and CpLs) low-pass configurations of two-element lossy L-Networks.

I do this by calculating the L-network L and C values required to match various loads to 50 ohms, given frequency and component Qs.  (Please refer to this post for the equations used for calculating L and C, and their derivation:
http://k6jca.blogspot.com/2018/09/l-networks-new-equations-for-better.html

The load impedances were calculated from gamma found while stepping around various Smith Chart "circles of constant SWR" in 2 degree increments.  L and C was then calculated for each impedance, and from these values (and their Qs), L-network power dissipation was calculated and plotted.

The code responsible for plotting the data is at the end of the listing.

 (Let me add -- unlike Dick, who is a true MATLAB expert, I am a MATLAB dilettante, which I am sure is reflected in my code, below.  MATLAB experts undoubtedly will find better ways of implementing what I have attempted to do).

% Plot power loss of a low-pass L-Network
% around circles of constant swr.
% k6jca

clear all
close all
clc

% The Givens:

F = 3.5e6;
Zo = 50;
Rsource = Zo;
Xsource = 0;
Q_L = 100;
Q_C = 100;

% Define the SWRs of the circles that will be plotted.
swr_array = [10 9 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5...
              3 2.5 2 1.8 1.6 1.4 1.2 1.1 1.05];

% Will plot a point every 2 degrees around the SWR,
% So that delta, in radians, is:
radian_delta = 2*pi/180;
w = 2*pi*F;


for swr_index = 1:length(swr_array)
    swr = swr_array(swr_index);
    rho = (swr-1)/(swr+1);
    for i = 1:181  % plotting 360 degrees
        % calculate the real and imaginary components of gamma.
        % (i.e. the x and y coordinates of each data point).
        real_gamma(swr_index,i) = rho*cos(radian_delta*(i-1));
        imag_gamma(swr_index,i) = rho*sin(radian_delta*(i-1));
        % from gamma, calculate Zload (and its real and imaginary
        % parts.
        Rld(i) = Zo*(1 - real_gamma(swr_index,i)^2 - imag_gamma(swr_index,i)^2)/((1 - real_gamma(swr_index,i))^2 + imag_gamma(swr_index,i)^2);
        Xld(i) = Zo*2*imag_gamma(swr_index,i)/((1 - real_gamma(swr_index,i))^2 + imag_gamma(swr_index,i)^2);
        Zld(i) = Rld(i) + 1i*Xld(i);

        % Given Zload, Zsource (which is 50 ohms here), Ql and Qc,
        % calculate the XL and XC values of the LsCp or CpLs network
        % than will match Zload to Zsource, using the equations
        % developed with MATLAB's Symbolic Match toobox (and
        % subsequently placed into the function calls, below).
        %
        % Note that for an LsCp network, these equations return
        % two XL and two XC terms.  Ditto for the CpLs network.
        [LsCp_XL(i,:), LsCp_XC(i,:)] = ser_L_shunt_C_Q_3(Rld(i),Xld(i),Rsource,Xsource,Q_L,Q_C);
        [CpLs_XL(i,:), CpLs_XC(i,:)] = shunt_C_ser_L_Q_3(Rld(i),Xld(i),Rsource,Xsource,Q_L,Q_C);

        % Of the four pairs of XL and XC terms returnd by the
        % two functions,(each function call returns two pairs),
        % only one pair is valid.  This pair will have
        % no imaginary components (that is, neither XL nor
        % XC is a complex number),
        % and both XL and XC will be positive.
        %
        % So let's first march through all 8 values of XL and XC and
        % identify any that are both real (not complex) AND positive.
        if(isreal(LsCp_XL(i,1)) && real(LsCp_XL(i,1)) > 0)
            itsValid(i,1) = 1;
        else
            itsValid(i,1) = 0;
        end
        if(isreal(LsCp_XL(i,2)) && real(LsCp_XL(i,2)) > 0)
            itsValid(i,2) = 1;
        else
            itsValid(i,2) = 0;
        end
        if(isreal(LsCp_XC(i,1)) && real(LsCp_XC(i,1)) > 0)
            itsValid(i,3) = 1;
        else
            itsValid(i,3) = 0;
        end
        if(isreal(LsCp_XC(i,2)) && real(LsCp_XC(i,2)) > 0)
            itsValid(i,4) = 1;
        else
            itsValid(i,4) = 0;
        end
        if(isreal(CpLs_XL(i,1)) && real(CpLs_XL(i,1)) > 0)
            itsValid(i,5) = 1;
        else
            itsValid(i,5) = 0;
        end
        if(isreal(CpLs_XL(i,2)) && real(CpLs_XL(i,2)) > 0)
            itsValid(i,6) = 1;
        else
            itsValid(i,6) = 0;
        end
        if(isreal(CpLs_XC(i,1)) && real(CpLs_XC(i,1)) > 0)
            itsValid(i,7) = 1;
        else
            itsValid(i,7) = 0;
        end
        if(isreal(CpLs_XC(i,2)) && real(CpLs_XC(i,2)) > 0)
            itsValid(i,8) = 1;
        else
            itsValid(i,8) = 0;
        end

        % Next, we are looking for "pairs" of XL and XC that are valid.
        % Only one pair should be valid.  The other three pairs should
        % be invalid.
        % For example, let's first check XL(1) and XC(1) for LsCp networks
        % and see if both are valid...
        LsCp(i) = 0;
        if(itsValid(i,1) && itsValid(i,3))    % XL(1), XC(1), LsCp netwrk
            LsCp(i) = 1;
            XL(i) = LsCp_XL(i,1);
            XC(i) = LsCp_XC(i,1);
        elseif(itsValid(i,2) && itsValid(i,4)) % XL(2), XC(2), LsCp netwrk
            LsCp(i) = 1;
            XL(i) = LsCp_XL(i,2);
            XC(i) = LsCp_XC(i,2);
        elseif(itsValid(i,5) && itsValid(i,7)) % XL(1), XC(1), CpLs netwrk
            XL(i) = CpLs_XL(i,1);
            XC(i) = CpLs_XC(i,1);
        elseif(itsValid(i,6) && itsValid(i,8)) % XL(2), XC(2), CpLs netwrk
            XL(i) = CpLs_XL(i,2);
            XC(i) = CpLs_XC(i,2);
        else
            % if there are no valid pairs...
            LsCp(i) = 2  % just a flag to id where the problem is.
            swr
            i
            disp('huh?');
        end

        % from XL and XC, calculate Ls and Cp
        Ls(i) = XL(i)/w;
        Cp(i) = 1/(XC(i)*w);

        % and calculate the S-parameters based upon whether the valid pair
        % is for an LsCp network or a CpLs network...
        if(LsCp(i))
            Spar = S_LsCp(Ls(i),Q_L,Cp(i),Q_C, Zo, F);
        else
            Spar = S_CpLs(Ls(i),Q_L,Cp(i),Q_C, Zo, F);
        end

        % With the s-parameters and Zload, calculate power loss
        Gp_Lnet(i) = powergain(Spar,Zo,Zld(i),'Gp');
        dBLoss(i) = 10*log10(Gp_Lnet(i));
        Percent_Loss(swr_index,i) = (1-Gp_Lnet(i))*100;

    end;

end

% Plot 3-d smith chart

% First, set Smith Chart parameters for smith_rab_v2.m,
% per its documentation:

LW=2;  % Line Width
SP.Rvalues=[0 10 25 50 100 250];   % SP.*= Smith Parameters
SP.Xvalues=[10 25 50 100 200 500];
SP.Zo=50;
SP.Nseg=60;
SP.LW  = 1;                 % LW= line width for Smith coordinates
SP.swr_circles = [];        % set to [] for no SWR circles
SP.LW_swr = 2;              % LW= line width for SWR Circles

SP.colors.grid   = [1 1 0]*0.5;       % coordinate lines
SP.colors.fill   = [0,0,0];           % background
SP.colors.inner_text = [1 1 0];       % inner labels
SP.colors.outer_text = [1 1 0];       % outer labels
SP.colors.swr    = [1 0 0];           % swr circle color
SP.colors.Q      = [0 0 1];           % Q contour color

SP.Q_contours    = [];
SP.Q_pts         = [];  % number of points in a Q contour
SP.LW_Q          = 1;

my_color_choice = [0.5 0.5 0.5];

% And now, PLOT!
handles.figure(4)     = figure('NumberTitle','Off',...
                        'Name', ' Match Space',...
                        'color',my_color_choice,...
                        'Position',[380 300 756 695]);
handles.Axes_Smith(4) = axes;
handles.smith(4)      = smith_rab_v2(handles.Axes_Smith(4),SP);
set(handles.Axes_Smith(4),'visible','on',...
                          'Xcolor',my_color_choice,...
                          'Ycolor',my_color_choice,...
                          'Zcolor',SP.colors.outer_text);
zlabel(['Gp Power Loss in PERCENT']);
title(handles.Axes_Smith(4),['Gp Percent Power Loss'],'Color',[1 1 0]);
colormap jet;

handles.cb(1)=colorbar;
set(handles.cb(1),'Color', SP.colors.outer_text,...
                  'FontWeight','Bold');
set(handles.cb(1).Label,'String','Percent Power Loss',...
                  'Color',SP.colors.outer_text,...
                  'FontWeight','Bold');

handles.surf(1,:) = surfc(real_gamma,imag_gamma,Percent_Loss);
set(handles.surf(1,1),'Edgecolor','none',...
                      'FaceAlpha',0.5,...
                      'LineStyle','none');
set(handles.surf(1,2),'Linewidth',2);  % nice fat contours lines!

sinfo = {['Ql : ',num2str(Q_L)],...
         ['Qc : ',num2str(Q_C)]};

ht=text(-1,1,15,sinfo);
set(ht,'Color',[0 1 0],'FontSize',12)
 

Published with MATLAB® R2018a

Plots from the code, above:

Downloading MATLAB functions for Plotting VNA Data:

Dick Benson's MATLAB functions and Utilities can be downloaded from his page in the MathWorks file exchange section:  Dick Benson MATLAB files.

Note that his enhanced Smith Chart function is in the "S-Parameter Utilities" (version 1.0.1) zip file.  A documentation PDF is included in the download.


Antenna Tuner Blog Posts:

A quick tutorial on Smith Chart basics:
http://k6jca.blogspot.com/2015/03/a-brief-tutorial-on-smith-charts.html

Plotting Smith Chart Data in 3-D:
http://k6jca.blogspot.com/2018/09/plotting-3-d-smith-charts-with-matlab.html

The L-network:
http://k6jca.blogspot.com/2015/03/notes-on-antenna-tuners-l-network-and.html

A correction to the usual L-network design constraints:
http://k6jca.blogspot.com/2015/04/revisiting-l-network-equations-and.html

Calculating L-Network values when the components are lossy:
http://k6jca.blogspot.com/2018/09/l-networks-new-equations-for-better.html

A look at highpass T-Networks:
http://k6jca.blogspot.com/2015/04/notes-on-antenna-tuners-t-network-part-1.html

More on the W8ZR EZ-Tuner:
http://k6jca.blogspot.com/2015/05/notes-on-antenna-tuners-more-on-w8zr-ez.html  (Note that this tuner is also discussed in the highpass T-Network post).

The Elecraft KAT-500:
http://k6jca.blogspot.com/2015/05/notes-on-antenna-tuners-elecraft-kat500.html

The Nye Viking MB-V-A tuner and the Rohde Coupler:
http://k6jca.blogspot.com/2015/05/notes-on-antenna-tuners-nye-viking-mb-v.html

The Drake MN-4 Tuner:
http://k6jca.blogspot.com/2018/08/notes-on-antenna-tuners-drake-mn-4.html

Measuring a Tuner's "Match-Space":
http://k6jca.blogspot.com/2018/08/notes-on-antenna-tuners-determining.html

Measuring Tuner Power Loss:
http://k6jca.blogspot.com/2018/08/additional-notes-on-measuring-antenna.html

Standard Caveat:

As always, I might have made a mistake in my equations, assumptions, drawings, or interpretations.  If you see anything you believe to be in error or if anything is confusing, please feel free to contact me or comment below.

And so I should add -- this information is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.