## Copyright (C) 2020 Reinhard Kotucha <reinhard.kotucha@web.de>
## 
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public License
## as published by the Free Software Foundation; either version 2
## of the License, or (at your option) any later version.
## 
## This program 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.  See the
## GNU General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
## MA 02110-1301, USA.

## -*- texinfo -*-
## @deftypefn {FunctionFile} {} smithchart()
## Plot smitchchart.
## @end deftypefn

# $Id: smithchart.m 1 2021-05-15 22:12:05Z reinhard $

function smithchart ()
  circlecolor = [.5 0 0];
  arccolor = [0.9 0.5 0.5];
  linewidth = 1.5;
  reso = 0.05; # Resolution.  Arcs are approximated by polygons.
  hold on;
  
  intersections = [
                   0.2  2;
                   0.5  2;
                   1    5;
                   2    5;
                   5   10;
                   10 500;
  ];

  for t = 1:rows(intersections)
    a = intersections(t,1); # arc
    b = intersections(t,2); # end point
    ## im(z) = 0
    plot([-1 1], [0,0])
    ## real
    z=a+[-b:reso:b].*i;
    plot((z-1)./(z+1), 'color', arccolor, 'linewidth', linewidth)
    ## imaginary
    z=[0:reso:b]+a*i;
    plot((z-1)./(z+1), 'color', arccolor, 'linewidth', linewidth)
    z=[0:reso:b]-a*i;
    plot((z-1)./(z+1), 'color', arccolor, 'linewidth', linewidth)
    labels(a)
  end
  plot(exp(([0:1:360].*i./180).*pi), ...
       'color', circlecolor, 'linewidth', linewidth)
  xlabel("Re Γ(z)");
  ylabel("Im Γ(z)");
  axis ([-1 1 -1 1], "square"); # aspect ratio
endfunction


function labels (a)
  if a<0
    label = sprintf("%.1f", a);
  else
    label = sprintf("%.d", a);
  end

  z=a+0i;
  G=(z-1)./(z+1);
  text(real(G), imag(G), ...
       label, "horizontalalignment", "right")

  z=0 + a*i;
  G=(z-1)./(z+1);
  text(real(G)+0.1*sin(G), imag(G)+0.05*cos(G), ...
       ["j" label], "horizontalalignment", "center")

  z=0 - a*i;
  G=(z-1)./(z+1);
  if imag(G) != -1 # Don't interfere with the Gamma axis
    text(real(G)+0.1*sin(G), imag(G)-0.05*cos(G), ...
         ["-j" label], "horizontalalignment", "center")
  endif
endfunction


# Local Variables:
#  mode: Octave
#  indent-tabs-mode: nil
#  coding: utf-8-unix
# End:
# vim:set tabstop=2 expandtab:
