function zviz(varargin);

%Modified by JBA to play nicer: Nov 20, 2015
% Upper case arguments are used to determine the syntax of the labels
% Concept: Extend the syntax to include output argument: w=s.^2
% The LHS arg is used for the labels.
%
% zviz has been modified so that the input line can contain Z=F(S), W=F(S), Z=F(W), etc 
% Assumed coordinates: W=u+iv; Z=x+iy; S=\sigma+j\omega
% S=F(Z) will label s=\sigma(x,y) +j\omega(x,y)
% W=F(S) label as u(\sigma,\omega), v(\sigma,\omega)
% F(S) will label s on the abcissa and ordinate
% F(W) will label as w=u+iv
%
% Example: W=S.^2
%


%ZVIZ  Visualize functions of a complex variable.
%   ZVIZ expr1 expr2 ... exprn or ZVIZ('expr1','expr2',...'exprn') plots
%   colored contour maps of the given expressions, which should involve
%   a complex variable, z.  Contour lines for the real and imaginary parts
%   are superimposed on a full color image using a color scheme developed
%   by John Richardson.  (http://home1.gte.net/jrsr/complex.html)
%   Red is real, blue is positive imaginary, green is negative imaginary,
%   black is small magnitude and white is large magnitude.  Branch cuts
%   appear as color discontinuities and coalescent contour lines.
%   Singularities appear as characteristic black and white patterns.
%
%   For example
%      zviz  z  2./z.^2  sqrt(z)  acosh(z)  atanh(z)  exp(-1./z.^2)
%   produces six figures with plots of the six expressions.
%
%   Regarding branch cuts, see pp. 79 and 86 of Abramowitz and Stegun
%   and Cleve's Corner, MATLAB News & Notes, Summer, 1998.
%
%                                   z = x + i*y
%   sqrt(z), log(z):                x <= 1       &  y == 0
%   asin(z), acos(z), atanh(z):     abs(x) >= 1  &  y == 0
%   atan(z), asinh(z):              x == 0       &  abs(y) >= 1
%   acsc(z), asec(z), acoth(z):     abs(x) <= 1  &  y == 0
%   acot(z), acsch(z):              x == 0       &  abs(y) <= 1
%   acosh(z):                       x <= 1       &  y == 0
%   asech(z):                  (x <= 0 | x >= 1) &  y == 0

set(0,'defaulttextfontsize',20);

%if nargin == 0
%vargargin{1}='z.^2';
%vargargin{2}= 'z';
%nargin = 2;
%end

warnings = warning;
warning('off')
for k = 1:nargin

   F = varargin{k};
   W=[];

   %look for = and strip it off; Return the output argument in W
   K=strfind(F,'=');
   if K ~= 0;
      W=F(1:K-1); F=F(K+1:end); %
	k=k+1; %skip to the next function on the arg list
   end
   [x,y] = meshgrid(-2:1/32:2);
   z = x + i*y;
   w = feval(inline(F),z);
 %  w = flipud(w); % srr fix y axis
%   w = fliplr(w); % jba fix x axis

%jba figure layering program (function defined below)
   layer(k);
   image(z2rgb(w))
   axis square
   tic = 1:32:129;
   xlabels = ['-2';'-1';' 0';' 1';' 2']; ylabels = [' 2';' 1';' 0';'-1';'-2'];
   set(gca,'xtick',tic,'xticklabel',xlabels,'ytick',tic,'yticklabel',ylabels);

   hold on
	if strfind(F,'Z') % W = F(x+jy) 
   Fs = regexprep(F,'Z','(x+jy)'); %if z) replace with s)
%   Fs = regexprep(Fs,'Z','z'); %if z) replace with s)
   title(['w = ',Fs],'fontsize',24)
   xlabel('x','fontsize',24)
   ylabel('jy','fontsize',24)
   
	elseif strfind(F,'S') % W = F(sigma+jomega)
%   Fs = regexprep(F,'S)','(\\sigma+j\\omega))'); %if z) replace with s)
%   Fs = regexprep(F,'S','(\\sigma+j\\omega)'); %if z) replace with s)
    Fs = regexprep(F,'S','s');
   title(['w = ',Fs],'fontsize',24)
   xlabel('\sigma','fontsize',24)
   ylabel('j\omega','fontsize',24)
   
	elseif strfind(F,'W')
   Fs = regexprep(F,'W','(u+jv)'); %if z) replace with s)
   title(['s = ',Fs],'fontsize',24)
   xlabel(['u'],'fontsize',24)
   ylabel(['iv'],'fontsize',24)
	end

   c = -4*pi:pi/4:4*pi; %line every pi/4 (45 degrees)
  contour(real(w),c,'k') %black=contour plots
  contour(imag(w),c,'w:') %white
%   contour(real(w),c) %color contours
%   contour(imag(w),c) %color contours

%plot a point
 %  plot(0.0,0.0,'wo');
   hold off

   drawnow
end

function [h]=layer(i);
% jont Allen for offset/layering figures for easy access on screen
h=figure(i);
P=get(h,'Position'); %left, bottom, width, height
U=get(h,'Units'); %system of units
%set(h,'Position',[001 524 560 420])
Q=[P(1)+30*(i-2),(i-1)*30,2*560,2*420];
set(h,'Position',Q);


function C = z2rgb(z);
%Z2RGB Complex variable color image.
%   Z2RGB(Z) maps a complex matrix into a full color image.
%   The mapping is due to John Richardson
%   See http://home1.gte.net/jrsr/complex.html.
%   Example:
%      [x,y] = meshgrid(-2:1/16:2); z = x+i*y;
%      image(z2rgb(F(z)))

r = abs(z);
a = sqrt(1/6)*real(z);
b = sqrt(1/2)*imag(z);
d = 1./(1+r.^2);
R = 1/2 + sqrt(2/3)*real(z).*d;
G = 1/2 - d.*(a-b);
B = 1/2 - d.*(a+b);
d = 1/2 - r.*d;
d(r<1) = -d(r<1);
C(:,:,1) = R + d;
C(:,:,2) = G + d;
C(:,:,3) = B + d;

