% CONVERGE1.M
%
% Let z be a complex number. We know z^n converges if ||z|| < 1, and that
% z^n goes to infinity if ||z|| > 1. However, if ||z|| = 1, then some
% interesting behavior occurs. We will examine this behavior graphically
% in the complex plane.
%
% Jessica Bernier        MATH215        Sep. 04, 2008

fprintf(['\nLet z be a complex number of the form z = a + bi,\n   where a and b are real numbers and i = sqrt(-1).\n']);
fprintf(['   Consider a and b in the range [-1,1] so that ||z|| = 1.\n']);

a = input('Enter the real component of the complex number (a): ');
while abs(a)>1
    a = input('Error: Please enter a value between -1 and 1: ');
end
b = sqrt(abs(1-a^2));
disp(['Given a = ', num2str(a), ', we need b = +/-', num2str(b), '.']);
disp(['The ratio (b^2)/(a^2) = ', num2str((b/a)^2), '.']);

nmax = 500; %This is the maximum value of n to use. Adjust nmax to show more or fewer points, if needed.
z1 = complex(a,b);
phi = angle(z1);
theta = phi/2/pi;
[Num,Denom] = rat(theta);
disp(['The initial angle (in polar coordinates) = 2pi * ', num2str(Num), '/', num2str(Denom), '.']);

figure;
plot(complex(cos(0:pi/1000:2*pi),sin(0:pi/1000:2*pi)),'k'); %Plot the unit circle in black, as a guide.
hold on; axis equal; grid on;
for n = 1 : nmax/5  %Plot the first 100 points in blue
    z = exp(i*phi*n); %The radius will always be 1.
    if isreal(z)==1 %If Matlab considers z to be real, it will graph in Cartesian coordinates. We need to manually declare z to be complex. This occurs for a = 1.
        plot(complex(z), '.b'); %complex(z) = z + 0i for real z.
    else
        plot(z, '.b');
    end
    pause(1/nmax) %Delay updating the figure so we can observe each new point being plotted.
end
hold on; axis equal; grid on;
plot(complex(cos(0:pi/1000:2*pi),sin(0:pi/1000:2*pi)),'k'); %Replot the unit circle in black, in case the figure was closed.
for n = nmax/5+1 : nmax*2/5  %Plot the next 100 points in red
    z = exp(i*phi*n); 
    if isreal(z)==1 
        plot(complex(z), '.r'); 
    else
        plot(z, '.r');
    end
    pause(1/nmax) 
end
hold on; axis equal; grid on;
plot(complex(cos(0:pi/1000:2*pi),sin(0:pi/1000:2*pi)),'k');
for n = nmax*2/5+1 : nmax*3/5  %Plot the next 100 points in green
    z = exp(i*phi*n); 
    if isreal(z)==1
        plot(complex(z), '.g');
    else
        plot(z, '.g');
    end
    pause(1/nmax)
end
hold on; axis equal; grid on;
plot(complex(cos(0:pi/1000:2*pi),sin(0:pi/1000:2*pi)),'k'); 
for n = nmax*3/5+1 : nmax*4/5  %Plot the next 100 points in magenta
    z = exp(i*phi*n); 
    if isreal(z)==1 
        plot(complex(z), '.m');
    else
        plot(z, '.m');
    end
    pause(1/nmax) 
end
hold on; axis equal; grid on;
plot(complex(cos(0:pi/1000:2*pi),sin(0:pi/1000:2*pi)),'k');
for n = nmax*4/5+1 : nmax  %Plot the next 100 points in cyan
    z = exp(i*phi*n);
    if isreal(z)==1
        plot(complex(z), '.c');
    else
        plot(z, '.c');
    end
    pause(1/nmax)
end
hold off
title('Behavior of z^n on the complex plane, for ||z|| = 1')
xlabel('Real component (a)')
ylabel('Imaginary component (b)')

disp('All done!');