function fcn(action)
%  set range= 1 to plot only z such that |f(z)|<1,
%            =-1 to plot only z such that Im f(z)>0
%            =-2 to plot only z such thatn Re f(z) > 0.
%            = 0 to plot all z
% jagged lines can be reduced by lowering the 
% resolution, but it takes longer.
% Warning: contour doesn't work well with discont fcn. Use NaN
% to remove pts where not cont. See e.g.s below
% To cut down on domain:e.g.:  w(imag(z) < 0)=NaN; 
% NEED: axis equal before imagesc and axis square after imagesc   

global Size
global xcenter
global ycenter
global resol
global range
global domain
global dimension 
global mmod
global modlines
global arglines
if nargin < 1
% Choices:
   domain=0; %plane=0, disk=1, ext(disk)=2, UHP=-1, RHP=-2
   range=0; % plane=0, disk=1, ext(disk)=2, UHP=-1, RHP=-2
   dimension=2; % 2d plot=2 3D plot =3
   Size=8; %width of plot (domain)
   resol=0.001; %relative grid size of domain
   mmod=1.8; % modulus lines in range: -mmod \le log|f| \le mmod
   modlines=20; % number of modulus lines, well -1.
   arglines=36; % number of argument lines
   xcenter=0; % x coordinate of center of plot
   ycenter=0 ; %y coordinate of center of plot
% end choices, except for the function.   
   fname = mfilename;
   callbackStr=[fname '(''zoomin'')'];
   uicontrol(...
        'Style','pushbutton',...
        'Units','normalized',...
        'Position', [.70 .02 .1 .03],...
        'Background','white',...
        'Foreground','black',...
        'String',' zoomin',...
        'tag','zoomin',...
        'Callback',callbackStr);
   callbackStr=[fname '(''zoomout'')'];
   uicontrol(...
        'Style','pushbutton',...
        'Units','normalized',...
        'Position', [.80 .02 .1 .03],...
        'Background','white',...
        'Foreground','black',...
        'String',' zoomout',...
        'tag','zoomout',...
        'Callback',callbackStr);
   callbackStr=[fname '(''quit'')'];
   uicontrol(...
        'Style','pushbutton',...
        'Units','normalized',...
        'Position', [.90 .02 .1 .03],...
        'Background','white',...
        'Foreground','black',...
        'String',' quit',...
        'tag','quit',...
        'Callback',callbackStr);
   callbackStr=[fname '(''2D/3D'')'];
   uicontrol(...
        'Style','pushbutton',...
        'Units','normalized',...
        'Position', [.60 .02 .1 .03],...
        'Background','white',...
        'Foreground','black',...
        'String',' 2D/3D',...
        'tag','2D/3D',...
        'Callback',callbackStr);
elseif strcmp(action,'zoomin')
   [xcenter,ycenter]=ginput(1);
   Size=Size/2;
elseif strcmp(action,'zoomout')
   Size=Size*2;
elseif strcmp(action,'quit')
   close;
   return
elseif strcmp(action,'2D/3D')
   dimension=5-dimension; 
end
   deltax=resol*Size;
   X = -Size/2.+xcenter:deltax:Size/2.+xcenter;
   Y=X-xcenter+ycenter;
   [x,y]=meshgrid(X,Y);
   z = x+1i*y;
% some functions to try
% w=z;
% w=(z-1i)./(z+1i); %Cayley Transform
  %w=(z-1).*(z+2).*(z-i);
%  w=(z-3).^2./(z.^2+4) ; 
%w=exp(z);
%    w=log(z);
%   w=log(z)+2*pi*1i;
%    w=z.^1.5
%   w=z.^4;
%     w=z.^(1/4);
%    w=exp(z.^3);

%    w=tan(z);
% w=(z+1./z)/2;
%    w=exp((z+1./z)/2);
%   w=.5*z.*(1+sqrt(1-1./z.^2));%
%     w=1i*(1i*z).^0.5-.5*1i;
%   w=exp(z+1./z);
%       w=(exp(z)+exp(-z))/2;
%      w=(w+1i)./(w-1i);    
%  w=(besselj(0,z).^2-besselj(1,z).^2).*z./sin(2*z);
%   w=(besselj(1,z).^2-besselj(3,z).^2).*z./sin(2*z);
%  w=1./z-1./sin(z);
%      w=besselj(0,z);
%Joukovski
%  w=z.*(1+sqrt(1-1./z.^2)); 
%   w=1./z;
%    w=z+1i*sqrt(1-z.^2);% add NaN statement below
%    w=z-1i*sqrt(1-z.^2);% add NaN statement below    
%     w=-z.^3;
%     w(imag(z) < 0)=NaN;
%   w(imag(z)==0 & real(z).^2 > 1)=NaN; %Kill where mod is not cont.
%  NOTE: if extraneous contours occur, decrease resol.   
%   Same map as sequence of compositions,
%      w=-1i*sqrt((1-z)./(1+z));w=(1+w)./(1-w);
%      w(imag(z)==0 & real(z).^2 > 1)=NaN; %Kill where mod is not cont.  
% to get arccos using J map:
%     w=-1i*log(w);
%    w=exp(15./z);
%    w=(z-1i/3)./(1+1i*z/3);
%    w=z.^(4);
%   w=sin(z)-z;
%  w=cos(z);
%    w=z./(z-(1+1i));
%     w=((z.^4-1)./(z.^4+1));
%   w=sqrt(z.^2+1);
%     w=z./(1-z);
%      w=(z.*sqrt(1+1./z.^2));
%      w=-1i*log(z+1i*sqrt(1-z.^2));      
w=exp((z+1)./(z-1));
%  M.C. Escher:
% the push forward of a polar grid
% by the map exp((1+i)z/2) is the same as the 
%pull back of the polar grid by the inverse map
%(1-i)log(w). The push forward of a rectangular grid
% is the same as the pull back of the polar grid by the map
% exp((1-i)logw)= w^{1-i).
%    w=exp((1-1i)*log(z)); 
% w=z.^(1-1i);w(imag(z)==0 & real(z) < 0 )=NaN;%Kill pts where mod is not continuous.
%   If we color the plane with a rectangular grid instead of polar, 
% then w would have the coloring of:
% w=exp(z);
%   w=sqrt(1i*z+.25)-0.5;%exterior of parabola
%    w=1i*w;
%    w=exp(z);%strip to half disk
%     w(abs(w)>1)=NaN;
%     w(imag(z)>pi)=NaN;
%     w(imag(z)<0)=NaN;
%       w=z;w(imag(z)<0)=NaN;
%Plot only points in the range:  
   if range==1
     w(abs(w)>1)=NaN; 
   elseif range==2
     w(abs(w)<1)=NaN;  
   elseif range==-1
     w(imag(w)<0)=NaN; 
   elseif range==-2
     w(real(w)<0)=NaN;
   end
   if domain==1
       w(abs(z)>1)=NaN;
   elseif domain==2
       w(abs(z)<1)=NaN;
   elseif domain==-1
       w(imag(z)<0)=NaN;
   elseif domain==-2
       w(real(z) < 0)=NaN;
   end    
   ww=angle(w)/pi;
   ww(imag(w)==0 & real(w) < 0)=-1; %fix red color at discontinuity
   axis equal
%Use the following instead of imagesc for angle lines.
%   contour(x,y,ww,-1:2/arglines:1);
% when angle=\pm pi all colors occur (contin approx of jump fcn)
  if dimension==2
    ww(isnan(ww))=2; %NaNs are white
    imagesc(X,Y,ww);  
    view(0,90)

    rotate3d off
  elseif dimension==3
    surf(x,y,log(abs(w)),ww); %NaNs are automatically white
    rotate3d on
    view(-37.5,30)
  end  
  shading interp
  set(gca,'YDir','normal');
  caxis([-1,2*modlines/arglines+1])%color=[-1,1] gray=[1,last]
  axis square
  colormap([hsv(arglines);gray(modlines)])
  hold on
  ww2=log(abs(w));
  ww3=(mmod-ww2)*modlines/(mmod*arglines)+1; %flip order so small=white
  if dimension==2
    contour(x,y,ww3,1:2/arglines:2*modlines/arglines+1,'LineWidth',2); 
    contour(x,y,abs(w),[1 1],'-k','LineWidth',2)
  elseif dimension==3
%    contour3(x,y,log(abs(w)),-1:2/arglines:2*modlines/arglines);  
    contour3(x,y,ww2,[0 0],'-k') 
    xlabel('real'), ylabel('imaginary'), zlabel('log-modulus')
  end

  
  hold off
% some other approaches:  
%   contour(x,y,t/pi,[1 1],'-r','LineWidth',2)
%   contour(x,y,angle(w)/pi,[-1 -1],'-r','LineWidth',2)
%polar version of domain (not range):
%   radius=4;
%   resol=0.03;
%   T=-pi:resol:pi;
%   R=0:resol:radius;
%   [r,t]=meshgrid(R,T);   
%   [x,y]=pol2cart(t,r);  
%   hc=colorbar;
%   get(hc);
%   set(hc,'Position',[.887 .1060 .03 .8167]);
%   set(hc,'YTickLabel',{'-pi','-pi/2','0pi','pi/2','pi, l2','l0','l-2'});
%   set(gca,'YDir','normal');hold off 
%   title('Color = arg f(z) and log|f(z)|')
%   xlabel('real'), ylabel('imaginary')
%   mp=max(max((abs(w))).'); 
%  if(isinf(mp))
%    mp=1.e+2;
%   end
%   mp=min(log(mp),2);
%   mn=min(min((abs(w))).');
%   if(isinf(mn));
%     mn=1.e-2;
%   end
%   mn=max(log(mn),-2);
%   scale=2/(mp-mn);
%   shift=(mn+mp)/(mp-mn);
% at some point fix this so that the unit circle is included 
% when appropriate.
%   ww2=haw.*(scale*(log(abs(w)+eps))-shift-2)+3;
%      colorbar('horizontal','XTickLabel',{'-pi','-pi/2','0','pi/2','-2','0','2'});
%      colorbar('vertical')
% To plot modulus repeating colors periodically:
%  problem: messy near log|w|=k pi
%    ww=haw.*(angle(exp(1i*(log(abs(w))+eps)))/pi-2)+2;
%    contour(x,y,ww,-1:.05:1)
%plot unit circle: problem with colors at k pi boundary
%  with the repeating colors approach.
end
