function out = dcf(fun,mag,semi,grid,z0)
%
% Visual display of "infinitesimal conformality": At each grid point
% the image of a (typically very small) square is drawn after magni-
% fication (mag  >> 1). The edges are color-coded to identify reflec-
% tions (e.g. 'conj(z)') and common rotations.
% If the function is (complex) differentiable, then in the limit (mag
% to infinity), the image of each square should again be a square,
% generally rotated and of different size, but w/ same orientation.
% (the derivative at each point is the complex number such that multi-
% plication by that number effects the corresponding "amplitwist", cf.
% Needham's book "Visual Complex Analysis".
%
% DCF(fun,mag,semi,grid,z0)
% DCF('z',1,3,10,0)
% DCF('conj(z)')
% DCF('z.^2')
% DCF('z.^3')
% DCF('sqrt(z.*conj(z)).*z',100)  
% DCF('(1+sqrt(z.*conj(z))/10).*z.^2',100)
% DCF('imag(z).*z',100)  
% DCF('1./z')
% DCF('log(z')
% DCF('exp(z')
% DCF('sqrt(z')
% DCF('(1+sqrt((z-3).*conj(z-3))/10).*(z-3).^2',100)
% DCF('sec(z)',1000,0.1,10,pi/2)
% DCF('i*z./(z.^2+1)',10)  %% needs more mag to see analyticity
% DCF                (self-demonstrating)
% 
% The roles and default values of the optional arguments are:
%
% fun      = 'z.^2' string that evaluates to a complex numeric when z complex
% mag      =    10  magnification (zooming) factor
% semidiam =     3  semidiameter of plot window
% grid     =    10  grid size 
% z0       = 0+i*0  coordinate(s) of center of window
%
%
% Author:  Matthias  Kawski. http://math.asu.edu/~kawski
%          All rights reserved.
% Date:    September 2002.
% Updates: None yet.

% /////////////////////////////////////
% Default values for optional arguments
% /////////////////////////////////////
%
if ( nargin < 1 )
   fun='z.^2'; 
end
polrat = isnumeric(fun);
% function given by formula w/ 'z' or as coeff's of poly / rational
if ( nargin < 2 )
    mag = 10; 
    %% even to avoid possible, common singularity at 0 which confuses
    %% the autoscaling feature (quadrilaterals should not overlap)
end
if ( nargin < 3 )
    semi = 3;
end
if ( nargin < 4 )
        grid = 10;
end
if ( nargin < 5)
       z0 = 0 + i * 0;
end

% ////////////////////////////////////////
% The main body of the color plot function
% ////////////////////////////////////////
%

close all 

dx    = 2*semi/grid;
x0    = real( z0 );
y0    = imag( z0 );
[x,y] = meshgrid( x0-semi+dx/2 : dx : x0+semi-dx/2,...
                  y0-semi+dx/2 : dx : y0+semi-dx/2 );            
z     = x+i*y;
oo    = ones(size(z));

% zooming factor, input
ddz   = dx / mag;

CTR   =  f(fun,z);
TR    =  f(fun,z+( 1+i)*ddz*ones(size(z))) - CTR;
TL    =  f(fun,z+(-1+i)*ddz*ones(size(z))) - CTR;
BR    =  f(fun,z+( 1-i)*ddz*ones(size(z))) - CTR;
BL    =  f(fun,z+(-1-i)*ddz*ones(size(z))) - CTR;

% take averages (diagonals) 
SE = (TR - BL) / 2;
SW = (TL - BR) / 2;
CR = abs((SW -i*SE)) < max( abs(SE),abs(SW) )/100;
bk = max(max(abs(TR)));
if bk < 1.01 *  min(min(abs(TR))); bk = 3 * bk; end;
Df = z2rgb( 3 / bk * TR / (i+i));

% rescale, zooming factor output
mag = dx / 2.0 / max(max(abs(0.00000000001+[TR,TL,BR,BL])));

TR    = z + mag * TR;
TL    = z + mag * TL;
BR    = z + mag * BR;
BL    = z + mag * BL;


p=figure;
set(p,'Name',['Derivative (?) of'])
set(p,'Units','normalized')
set(p,'Position',[0,0.03,1,0.9])
hold on;
view(2)
axis equal
axis([x0-semi,x0+semi,y0-semi,y0+semi])
% % 	plot(real(TR(:)),imag(TR(:)),'or');
% % 	plot(real(BR(:)),imag(BR(:)),'*g');
for i=1:size(z,1)
    for j=1:size(z,2)
        if CR(i,j)
            fill(real([TR(i,j),TL(i,j),BL(i,j),BR(i,j),TR(i,j)]),...
                 imag([TR(i,j),TL(i,j),BL(i,j),BR(i,j),TR(i,j)]),...
                 Df(i,j,:));
        end
        
% %             plot(real([TR(i,j),TL(i,j),BL(i,j),BR(i,j),TR(i,j)]),...
% %                  imag([TR(i,j),TL(i,j),BL(i,j),BR(i,j),TR(i,j)]),'b-');
        p1=plot(real([TR(i,j),TL(i,j)]),...
             imag([TR(i,j),TL(i,j)]),'g-');
        p2=plot(real([TL(i,j),BL(i,j)]),...
             imag([TL(i,j),BL(i,j)]),'b-');
        p3=plot(real([BL(i,j),BR(i,j)]),...
             imag([BL(i,j),BR(i,j)]),'k-');
        p4=plot(real([BR(i,j),TR(i,j)]),...
             imag([BR(i,j),TR(i,j)]),'r-');
        set(p1,'LineWidth',2.5)
        set(p2,'LineWidth',2.5)
        set(p3,'LineWidth',2.5)
        set(p4,'LineWidth',2.5)
     end;
 end;

% % 	if polrat
% %         if size(fun,1) == 1
% %             w=polyval(fun,z);
% %         else
% %             w=polyval(fun(1,:),z)./polyval(fun(2,:),z);
% %             if ratdiff
% %                 w=w-polyval(fun2(1,:),z)./polyval(fun2(2,:),z);
% %             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.
% % 	
% % 	p=surf(x,y,zeros(size(x)));
% % 	set(p,'CData',xy2rgb(u,v,color,brite),...
% %           'LineStyle','none',...
% %           'FaceColor','interp')
% % 	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
out=p;
return

% 	p=surf(x,y,zeros(size(x)));
% 	set(p,'CData',xy2rgb(u,v,color,brite),...
%           'LineStyle','none',...
%           'FaceColor','interp')

%--------------------------
% 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=z2rgb(z)
color = 1;
brite = 0;
% assume hue        h between 0 and 6
% assume brightness b between 0 and 3
eps = 0.00000001;
z = 2 * z / (-sqrt(3)+i);
r   = abs(z)+eps;
up  = (1-tanh(log(r)/color-brite))/2;
h1  = g(up).*real(z)./r/0.9; % the /0.9 makes the max(h2,h2)=unit length
h2  = g(up).*imag(z)./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