%
%  code02Wa.m -- djm -- 13 sept 2006
%
%  complex function mapping demo
%

%  parameters
L = 4;  N = 24;  

%  menu
choice = [];  
while (isempty(choice)==1)
  disp(' ')
  disp('  map demo:')
  disp('    0) initialize figure')
  disp('    1) line')
  disp('    2) box')
  disp('    3) circle')
  disp(' ')

  choice = input(' choice:  ');
end

%  default figure #
if(exist('figN')==0);  figN = 1;  end
    
switch choice
    %  initialize figure
    case {0}
        figN = input(' figure #:  ');
        figure(figN+10);  clf;  hold on
        axis equal;  axis([-1 1 -1 1]*4)
        figure(figN   );  clf;  hold on
        axis equal;  axis([-1 1 -1 1]*4)
    %  line map
    case {1}
        figure(figN)
        disp('');
        disp(' click to choose midpoint & endpoint');
        disp('');
        
        z0 = [];  z1 = [];
        while (isempty(z0))
            in = floor(ginput(1)/(pi/25))*(pi/25);

            z0 = in(1) + i*in(2);
            plot(in(1),in(2),'k*')

            in = floor(ginput(1)/(pi/25))*(pi/25);

            z1 = in(1) + i*in(2);
            a0 = z1-z0;
            plot(in(1),in(2),'ko')

            disp(['  centre = ' num2str(z0) ';  slope = ' num2str(a0)]);
        end
        
        % plot pre-image
        s = -1:2/N:1;
        pts = z0 + a0.*s;
        plot(real(pts),imag(pts),'r')
    %  box map
    case {2}
        figure(figN)
        disp('');
        disp(' click to choose midpoint & corner');
        disp('');
        
        z0 = [];  z1 = [];
        while (isempty(z0))
            in = floor(ginput(1)/(pi/25))*(pi/25);

            z0 = in(1) + i*in(2);
            plot(in(1),in(2),'k*')

            in = floor(ginput(1)/(pi/25))*(pi/25);

            z1 = in(1) + i*in(2);
            a0 = z1-z0;
            plot(in(1),in(2),'ko')

            disp(['  centre = ' num2str(z0) ';  diag = ' num2str(a0)]);
        end
        
        % plot pre-image
        s = -1:2/N:1;
        pts = [z0 +     real(a0) - i*imag(a0).*s , ...
                 z0 - i*imag(a0) -   real(a0).*s , ...
                 z0 -   real(a0) + i*imag(a0).*s , ...
                 z0 + i*imag(a0) +   real(a0).*s];
    %  circle map
    case {3}
        figure(figN)
        disp('');
        disp(' click to choose centre & circle point');
        disp('');
        
        z0 = [];  z1 = [];
        while (isempty(z0))
            in = floor(ginput(1)/(pi/25))*(pi/25);

            z0 = in(1) + i*in(2);
            plot(in(1),in(2),'k*')

            in = floor(ginput(1)/(pi/25))*(pi/25);

            z1 = in(1) + i*in(2);
            a0 = z1-z0;
            plot(in(1),in(2),'ko')

            disp(['  centre = ' num2str(z0) ';  radius = ' num2str(abs(a0))]);
        end
        
        % plot pre-image
        s = -1:2/N/4:1;
        pts = [z0 + abs(a0)*exp(i*2*pi*s)];
end

switch choice
    case {1,2,3}

        %  plot pre-image (pts = vector of complex pts)
        plot(real(pts),imag(pts),'r')
        title('\bf pre-image, z plane')
        xlabel('\bf Re-axis')
        ylabel('\bf Im-axis')

        %  TYPE map here
        f = pts.^3 - 2*pts.^2 - 2*pts.^1 - 3;
        
        %  plot image
        figure(figN+10)
        plot(real(f(1)),imag(f(1)),'ko')
        plot(real(f),imag(f),'b')
        psize = max(abs(f));
        axis equal;  axis([-1 1 -1 1]*1.1*psize);
        title('\bf image, w=f(z) plane')
        xlabel('\bf Re-axis')
        ylabel('\bf Im-axis')
end