function zviz(varargin);
% version 3.0 Jun 4, 2019: Fix S=log(W)
% version 2.0 Oct 9, 2016
% Modified by JBA to fix x and y labels and to not crash when printing with Octave (Oct 8, 2016)
% 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:
% The RHS arg is used for the labels: |w|(s) maps to intensity; angle(w)(s) maps to color
% w=s-sqrt(j) has zero (dark spot) at sqrt(j)=e^(j\pi/2) = (1+j)/sqrt(2)
% Input syntax:
% zviz has been modified so that the input line can contain Z=F(S), W=F(Z), Z=F(S), etc 
% Assumed coordinates: W=u+iv; Z=x+iy; S=\sigma+j\omega
% S=F(Z) will label intensity: |F|(x,y) and hue: angle(F)(x,y) 
% W=F(S)  will label intensity: |F|(sigma,omega) and hue: angle(F)(sigma,omega) 
%
% Example: W=Z-sqrt(j) will have zero (dark spot) at (1+j)/sqrt(2)
% Principle value of phase (hue) is at \pm pi


%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,'defaultAxesFontSize',14)
set(0,'defaulttextfontsize',14);
 set(0,'defaulttextinterpreter','none')
% set(0,'defaulttextinterpreter','latex')

nargin1=nargin;

if nargin == 0
varargin{2}='Z.^2';
varargin{1}= 'Z';
nargin1 = 2;
end

warning('off'); %default off

for k = 1:nargin1

   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']; % SRR 10/9/16: Labels are right
   set(gca,'xtick',tic,'xticklabel',xlabels,'ytick',tic,'yticklabel',ylabels);
   set(gca,'fontsize',16); % SRR 10/9/16 make tick labels larger

   hold on
	if strfind(F,'Z') % S = F(x+jy) 
   Fs = regexprep(F,'Z','z'); %F contains the function; if function(Z) replace with 'z') %JBA Jun4, 2019
   title(['s = ',Fs],'fontsize',24)   %JBA Jun4, 2019
   %SRR 10/9/16: Each element of the image plot is w(z) calculated from the 
   %matrix z, which is arranged like the z plane (e.g. top-left is 
   %z = -2+2i, top-right is z = 2+2i, bottom left is z = -2-2i and
   %bottom-right is z=2-2i). This can be checked in the command window.
   %Therefore, the proper labels are x and iy (the domain). The matrix w(z)
   %that we plot using the command "image" has the following elements: 
   %top-left is w(-2+2i), top-right is w(2+2i), bottom left is w(-2-2i) 
   %and bottom-right w(2-2i) 
   xlabel('x','fontsize',24)
   ylabel('jy','fontsize',24)
   
%%% Change title
	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)
    if strfind(W,'Z'),
    W = regexprep(W,'Z','R+Xj'); %if z) replace with s)
    title([W,' = ',F],'fontsize',24)  %no latex conversion

    else 
    Fs = regexprep(F,'S','s');
   title(['w = ',Fs],'fontsize',24)  %JBA Jun4, 2019
    end
   
   xlabel('\sigma','fontsize',24)
   ylabel('j\omega','fontsize',24)
   
	elseif strfind(F,'W') %ck for function of W, -> (u+vj)
   Fs = regexprep(F,'W','(u+jv)'); %if z) replace with s)
   title(['s = ',Fs],'fontsize',24)
   xlabel(['u'],'fontsize',24)
   ylabel(['jv'],'fontsize',24)
   
%replace "Z =" with "R+Xj ="
   Fs = regexprep(W,'Z','R+Xj'); %if z) replace with s)
   title([Fs],'fontsize',24)
   xlabel(['\sigma'],'fontsize',24)
   ylabel(['\omega j'],'fontsize',24)

	end

%remove double spaces in title
% not working, not sure why
%   Fs = regexprep(F,'  ',''); %if z) replace with s)
%   title([Fs],'fontsize',24)

   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;
grid on;

%Prepare to make eps figure
%   drawnow
   h = gcf; 
   % SRR 10/9/16: if you use gcf instead of handle, when you turn the 
   %second figure off, the first figure becomes gcf
%set(h, 'visible', 'off'); %work-around for printing with octave
 drawnow; pause(.3); %This delay seems to be critical for proper figures
 print -depsc figure
%set(h, 'visible', 'on'); %work-around for printing with octave
end

function [h]=layer(k) %SRR 10/9/16: swapped i for k (i is imaginary number)
% jont Allen for offset/layering figures for easy access on screen
%SRR 10/9/16: Depending on your screen size, the direct use of pixel
%numbers can actually make the figures harder to access (FYI)
h=figure(); % SRR 10/9/16: figure() means we initialize new fig without overwriting currently open figs
P=get(h,'Position'); %left, bottom, width, height
U=get(h,'Units'); %system of units
%set(h,'Position',[001 524 560 420])
scl=1.3; %scale size of figure on screen (pixels)
Q=[P(1)+30*(k-2),(k-1)*30,scl*560,scl*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;

