function out = ccf(fun,fun2,semi,grid,z0,col,z0f1,z0f2)
%              ccf(fun,fun2,semi,grid,z0,col,z0f1,z0f2)
% or           ccf(fun,     semi,grid,z0,col,z0f1)
%
% modified in argentina! Sept 17, 2002
% updated April 2006
%
% CCF(fun,semi,grid,z0,col)
%
% Domain coloring of functions from complex plane to itself. Allows for
% symbolic function definitions (use quotes, careful with ./ and .^ etc.),
% polynomials (eneterd as row vector), rational funcitons (martix with
% two rows of same length -- pad w/ zeros as needed!), and for the
% difference of two functions, used e.g. for the error term in Taylor
% and Laurent expansions. MOst dramatic are animations with e.g. 
% series approximations of increasingly higher order, or root-locus
% diagrams.
% Main known bugs are due to 2002 addition to allow for difference 
% of two functions (e.g. one symbolic, and other numeric polynomial 
% = row vector) ... the program is supposed to pick this up automati-
% cally, but not all shifts of subsequent arguments seem to be correct,
% likely this appears only in pre-2002 examples that were not updated???
%
% CCF('z.^2',3,111,0)
% CCF('exp(1./z)')
% CCF('z.*sin(1./z)',.3,222,0,[3,0])
% CCF([1 0 0 1])
% CCF([1 0 0 1;1 0 0 -1])
% CCF                (self-demonstrating)
% 
% New experiments:
% Successive (symmetric) Laurent expansion
% N=40;r=2;R=3;C=[[-(1/R).^[N+1:-1:1],r.^[0:N-1]];[zeros(1,N),1,zeros(1,N)]];
%   for k=0:N; ccf(C(:,N+1-k:N+1+k),2*abs(R)); pause(0.3); end;
% Remainder terms for the same:
% for k=0:N; ccf(C(:,N+1-k:N+1+k),[0,1,-r-R;1,-r-R,r*R],2*abs(R));pause(0.3);end;
% 
% Taylor approximations of transcendental functions
% for N=[1:30]; ccf([fliplr(1./cumprod(1:N)),1],10,111,-3); pause; end;
%% WATCH CAREFULLY: some zeros jettison off the semi-disk to the right,
%% MAA Monthly Dec 2005, p.894 -- zeros asymptotically to left of parabola!
% for N=[1:30]; ccf([fliplr(real(1./cumprod(i*[1:N]))),1],10,111,-3); pause; end;
% for N=[1:30]; ccf([fliplr(imag(1./cumprod(i*[1:N]))),0],10,111,-3); pause; end;
% nice mistake -- lots of double roots for "obvious reasons"
% for N=[1:30]; ccf([fliplr(imag(1./cumprod(i*[1:N]))),1],10,111,-3); pause; end;
%
% increasing control gain k of transfer functions... (need better examples)
% for k=[-10:0.2:10] ; ccf([1 0 9;1 k 9],10); pause(0.01); end
%
% For non-MATLAB users: Exert caution with proper syntax
% for array operations:  .*  , ./  , and .^
% 
% The second thru fourth argument are optional. 
% Their roles and default values are:
%
% semidiam =     3  semidiameter of plot window
% grid     =   256  grid size 
% z0       = 0+i*0  coordinate(s) of center of window
% col      = [1,0]  colors and briteness:
%          if col(1) > 1 the range of colored magnitudes,
%                  i.e. between white and black is wider.
%          if col(2) > 0 then white extends to larger
%                  magnitudes, while if col(2)<0 then
%                  black apears for smaller magnitudes
%
% More basic examples: 
% z, conj(z), i*z, 1/z, real(z), z.^2, z.^3
% 1./(z.^2-1), (z.^2+1)./(z.^2-1)
% finite Taylor and Laurent approximations
% exp(z), sin(z), cos(z)
% exp(1/z), sin(1/z)
% sqrt(z), -sqrt(z), log(z), sqrt(z^2-1)
%
% Author:  Matthias  Kawski. http://math.la.asu.edu/~kawski
%          All rights reserved.
% Date:    November 23, 1999.
% Updates: April 2006: redefined thisfig to show succ approx in same window
%          Sept 2002: Expansion point of numeric polynomials added
%          July 2002: Extended functionality to handle polynomials specified
%                     as usual as vectors of complex numerical coefficients,
%                     and rational functions entered as 2 x n matrix. 

% /////////////////////////////////////
% Default values for optional arguments
% /////////////////////////////////////
%
if ( nargin < 1 )
   fun='sin(z.^3)'; 
end
polrat = isnumeric(fun);
% function given by formula w/ 'z' or as coeff's of poly / rational
if ( nargin < 2 )
    semi = 3;
    ratdiff = 0;
else
    ratdiff = isnumeric(fun) & (prod(size(fun2)) >1 );
end
% difference of two rational functions shifts all later args
if ratdiff 
    if ( nargin < 3 )
        semi = 3;
    end;
   if ( nargin < 4 )
        grid = 111;
   end
   if ( nargin < 5)
       z0 = 0 + i * 0;
   end
   if ( nargin < 6)  
      color = 1;
      brite = 0;
   else
      color = col(1);
      brite = col(2);      
   end
   if ( nargin < 7)
       z0f1 = 0;
   end
   if ( nargin < 8)
       z0f2 = 0;
   end
else
    % shift the names of all later arguments. REVERSE porder critical.
   if ( nargin < 6)
       z0f1 = 0;
   else
       z0f1= col;
   end
   end
   if ( nargin < 5)  
      color = 1;
      brite = 0;
   else
      color = z0(1); 
      brite = z0(2);      
   end
   if ( nargin < 4)
       z0 = 0 + i * 0;
   else
       z0 = grid;
   end
   if ( nargin < 3 )
        grid = 256;
    else
        grid = semi;
   end
   if ( nargin < 2 )
       semi = 3;
   else 
       semi = fun2;
   end;
   
% ////////////////////////////////////////
% The main body of the color plot function
% ////////////////////////////////////////
%
dx    = semi*2/grid;
x0    = real( z0 );
y0    = imag( z0 );
[x,y] = meshgrid( x0-semi : dx : x0+semi,...
                  y0-semi : dx : y0+semi );            
z     = x+i*y;
oo    = ones(size(z));

if polrat
    if size(fun,1) == 1
        w=polyval(fun,z-z0f1);
    else
        w=polyval(fun(1,:),z-z0f1)./polyval(fun(2,:),z-z0f1);
        if ratdiff
            w=w-polyval(fun2(1,:),z-z0f2)./polyval(fun2(2,:),z-z0f1);
        end;
    end
else
    w     = f(fun,z);
end
u     = real(w);
v     = imag(w);

cut = 100;
u   = max(min(u,cut*oo),-cut*oo); % to prevent overflow, color
v   = max(min(v,cut*oo),-cut*oo); % will be black anyhow.

thisfig = gcf;
set(thisfig,'Units','normalized')
set(thisfig,'Position',[0,0.03,1,0.9])

p=surf(x,y,zeros(size(x)));
set(p,'CData',xy2rgb(u,v,color,brite),...
      'LineStyle','none',...
      'FaceColor','interp')
view(2)
axis equal
axis([x0-semi,x0+semi,y0-semi,y0+semi])
if isnumeric(fun); 
    if size(fun,1)==1 
        title(strcat('pol = [',num2str(fun),']'))
    elseif ~ratdiff
        title(strcat('rational = [',num2str(fun(1,:)),'] / [',num2str(fun(2,:)),']'))
    else title('difference of rationals')
    end
else title(fun) 
end
% figure(thisfig)  %% leave this for the calling procedure...
out = thisfig;
return

%--------------------------
% Subfunctions
%--------------------------
function freturn = f(fun,z)
freturn = eval(fun);
return    
%
function greturn = g(x)
% input and output both in [0..1]. Almost a sawtooth.
% g(0)=g(1)=0, and g(1/2)=1 (almost)
greturn = 1.01*ones(size(x))-...
   sqrt((2*x-ones(size(x))).^2+0.01*ones(size(x)));
return    
%
function out=xy2rgb(u,v,color,brite)
% assume hue        h between 0 and 6
% assume brightness b between 0 and 3
eps = 0.000001;
r   = sqrt(u.^2+v.^2+eps);
up  = (1-tanh(log(r)/color-brite))/2;
h1  = g(up).*u./r/0.9; % the /0.9 makes the max(h2,h2)=unit length
h2  = g(up).*v./r/0.9; % the /0.9 makes the max(h1,h2)=unit length
%
q=5; % if q is too small the cone does not fit inside the RGB-cube
out(:,:,1)=min(1,max(0,(-sqrt(3)*h1+h2)/q+up));
out(:,:,2)=min(1,max(0,( sqrt(3)*h1+h2)/q+up));
out(:,:,3)=min(1,max(0,          -2*h2 /q+up));
return