%% Example Of A Mobius Map: w = -(Z-A)/(Az-1), -1<A<1
%% CSUMS Summer 2010 A.O. Hausknecht 
## Modifed 9/26/2010 to display A
%% Written using OCTAVE
function plotZMobius_CircularGrid(A = 0.5)
printf("A = %g\n", A);
Astr = sprintf("A = %g", A);

clf;

%% Create A Plot Of The Circular Grid
subplot (1, 2, 1)
title ("Grid Of Truncated Sectors")
hold on

%% Draw filled of truncated sectors
plotFilledSector(0*pi/4, 1, 1*pi/4, 1.25, [0.8 0.8 0.8])
plotFilledSector(1*pi/4, 1, 2*pi/4, 1.25, [0.8 0 0])

plotFilledSector(0*pi/4, 0.5, 1*pi/4, 1, [0.9 0.9 0.9])
plotFilledSector(1*pi/4, 0.5, 2*pi/4, 1, [1 0 0])
plotFilledSector(2*pi/4, 0.5, 3*pi/4, 1, [0 1 0])
plotFilledSector(3*pi/4, 0.5, 4*pi/4, 1, [0 0 1])
plotFilledSector(4*pi/4, 0.5, 5*pi/4, 1, [1 1 0])
plotFilledSector(5*pi/4, 0.5, 6*pi/4, 1, [0 1 1])
plotFilledSector(6*pi/4, 0.5, 7*pi/4, 1, [1 0 1])
plotFilledSector(7*pi/4, 0.5, 8*pi/4, 1, [0.5 0.5 1])

plotFilledSector(0*pi/4, 0.25, 1*pi/4, 0.5, [0.7 0.7 0.7])
plotFilledSector(1*pi/4, 0.25, 2*pi/4, 0.5, [0.8 0 0])
plotFilledSector(2*pi/4, 0.25, 3*pi/4, 0.5, [0 0.8 0])
plotFilledSector(3*pi/4, 0.25, 4*pi/4, 0.5, [0 0 0.8])
plotFilledSector(4*pi/4, 0.25, 5*pi/4, 0.5, [0.8 0.8 0])
plotFilledSector(5*pi/4, 0.25, 6*pi/4, 0.5, [0 0.8 0.8])
plotFilledSector(6*pi/4, 0.25, 7*pi/4, 0.5, [0.8 0 0.8])
plotFilledSector(7*pi/4, 0.25, 8*pi/4, 0.5, [0.3 0.3 0.8])


%% Plot the grid's circles 
t = linspace(0,2*pi, 50);
for r =  [0.25 2]
  plot( r*cos(t),  r*sin(t),'b');
end
plot( cos(t),  sin(t),'k', 'linewidth', 3);

%% Plot the grid's radial lines
r = 0:.1:2;
for theta = 0: pi/4.0: 2*pi
  plot( r*cos(theta), r*sin(theta), 'r')
end
%% Set the axes domain, range, and scaling.
axis('equal')
xyMax = 1.25;
axis([-xyMax xyMax -xyMax xyMax])

%% Plot Unit Tangents
plotUnitTangents(1 ,  'magenta', 'yellow')
plotUnitTangents((1+i)/sqrt(2) , 'green', 'cyan')

hold off

%% Create A Plot The Image Of The Circular Grid
subplot (1, 2, 2)
title (["Image Of Grid Under\nw = -(z-A)/(Az-1), ", Astr])
hold on

%% Draw images of filled truncated sectors
plotImageOfFilledSector(0*pi/4, 1, 1*pi/4, 1.25, A, [0.8 0.8 0.8])
plotImageOfFilledSector(1*pi/4, 1, 2*pi/4, 1.25, A, [0.8 0 0])

plotImageOfFilledSector(0*pi/4, 0.5, 1*pi/4, 1, A, [0.9 0.9 0.9])
plotImageOfFilledSector(1*pi/4, 0.5, 2*pi/4, 1, A, [1 0 0])
plotImageOfFilledSector(2*pi/4, 0.5, 3*pi/4, 1, A, [0 1 0])
plotImageOfFilledSector(3*pi/4, 0.5, 4*pi/4, 1, A, [0 0 1])
plotImageOfFilledSector(4*pi/4, 0.5, 5*pi/4, 1, A, [1 1 0])
plotImageOfFilledSector(5*pi/4, 0.5, 6*pi/4, 1, A, [0 1 1])
plotImageOfFilledSector(6*pi/4, 0.5, 7*pi/4, 1, A, [1 0 1])
plotImageOfFilledSector(7*pi/4, 0.5, 8*pi/4, 1, A, [0.5 0.5 1])

plotImageOfFilledSector(0*pi/4, 0.25, 1*pi/4, 0.5, A, [0.7 0.7 0.7])
plotImageOfFilledSector(1*pi/4, 0.25, 2*pi/4, 0.5, A, [0.8 0 0])
plotImageOfFilledSector(2*pi/4, 0.25, 3*pi/4, 0.5, A, [0 0.8 0])
plotImageOfFilledSector(3*pi/4, 0.25, 4*pi/4, 0.5, A, [0 0 0.8])
plotImageOfFilledSector(4*pi/4, 0.25, 5*pi/4, 0.5, A, [0.8 0.8 0])
plotImageOfFilledSector(5*pi/4, 0.25, 6*pi/4, 0.5, A, [0 0.8 0.8])
plotImageOfFilledSector(6*pi/4, 0.25, 7*pi/4, 0.5, A, [0.8 0 0.8])
plotImageOfFilledSector(7*pi/4, 0.25, 8*pi/4, 0.5, A, [0.3 0.3 0.8])


%% Plot images of the circles under w = -(z-A)/(Az-1)
t = linspace(0,2*pi, 100);
for r =  [0 2]
  z = r*(cos(t)+sin(t)*i);
  w = mobius(z,A);
  plot(real(w), imag(w),'b');
end
z = (cos(t)+sin(t)*i);
w =  mobius(z,A);
plot( real(w), imag(w),'k', 'linewidth', 3);

%% Plot images of the radial lines under w = mobius(z,A)
r = 0: .1: 2; 
for theta = 0: pi/4.0: 2*pi
  z = r*(cos(theta)+sin(theta)*i);
  w =  mobius(z,A);
  plot(real(w), imag(w), 'r')
end
%% Set the axes domain, range, and scaling.
axis('equal')
axis([-xyMax xyMax -xyMax xyMax])

%% Plot the images of the Unit Tangents
plotImagesOfUnitTangents(1 ,  A, 'magenta', 'yellow')
plotImagesOfUnitTangents((1+i)/sqrt(2) , A, 'green', 'cyan')

hold off
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Helper functions used above
function w = mobius(z, A)
         w = -(z-A)./(A*z-1);
end 
  
function plotFilledSector(t0, r0, t1, r1, color)
%% Fill a truncated sector
%% 1. Generate points along the bottom edge 
%%    which is part of a circle.
t = linspace(t0, t1, 10);
x1 = r0*cos(t); y1 = r0*sin(t);
%%
%% 2. Generate points along the top edge 
%%    in revese order and generate the points
t = fliplr (t);
x2 = r1*cos(t); y2 = r1*sin(t);
%%
%% 3. Now collect all the points of the regions
%%    boundary together in counter-clockwise order!
x = cat(2, x1, x2); y = cat(2, y1,y2);
%%
%% 4. Finally, fill the region
fill (x, y, color)
end

function plotImageOfFilledSector(t0, r0, t1, r1, A, color)
%% Fill the image of truncated sector
%% 1. Generate points along the bottom edge 
%%    which is part of a circle.
t = linspace(t0, t1, 20);
z = r0*( cos(t)+sin(t)*i);
w = mobius(z,A);
x1 = real(w); y1 = imag(w);
%%
%% 2. Generate points along the right edge 
%%    which is part of a circle.
r =   linspace(r0, r1, 20);
z = r*( cos(t1)+sin(t1)*i);
w = mobius(z,A);
x2 = real(w); y2 = imag(w);
%%
%% 3. Generate points along the top edge 
%%    in revese order.
t = fliplr (t);
z = r1*( cos(t)+sin(t)*i);
w = mobius(z,A);
x3 = real(w); y3 = imag(w);
%%
%% 4. Generate points along the left edge 
%%    in revese order.
r = fliplr (r);
z = r*( cos(t0)+sin(t0)*i);
w = mobius(z,A);
x4 = real(w); y4 = imag(w);
%%
%% 3. Now collect all the points of the region's
%%    boundary together in counter-clockwise order!
x = cat(2, x1, x2, x3, x4); y = cat(2, y1, y2, y3, y4);
%%
%% 4. Finally, fill the region with green
fill (x, y, color)
end

function plotUnitTangents(z, color1, color2)
%% Draw Unit Tangent vectors at z
%%
%% 1. Draw the radial tangent to at z.
dz = z/abs(z);
a = [real(z) real(z + dz)]; b = [imag(z) imag(z + dz)];
line(a, b, 'linewidth',3, 'color', color1);
%%
%% 2. Draw the tangent dz to the circle at z.
%%    Note: Use J(x+yi) = -y+xi to obtain perpindicular vector
dz = -imag(dz) + real(dz)*i;
a = [real(z) real(z + dz)]; b = [imag(z) imag(z + dz)];
line(a,b, 'linewidth',3, 'color', color2)
end

function plotImagesOfUnitTangents(z, A, color1, color2)
%% Plot images of the Unit Tangent vectors at w = -(z-A)/(Az-1)
w = -(z-A)/(A*z-1);
%% 1. Use the derivative to plot the image of the radial tangent  
%%    to the circle:  Note: dw =(A^2 - 1)/(Az-1)*dz
dz = z/abs(z); dw =(A^2 - 1)/(A*z-1)^2*dz;
a = [real(w)  real(w+dw)]; b = [imag(w) imag(w+dw)];
line(a,b, 'linewidth', 3, 'color', color1);
%%
%% 2. Draw the tangent dz to the circle at z.
%%    Note: Use J(x+yi) = -y+xi to obtain perpindicular vector
dz = -imag(dz) + real(dz)*i; dw =(A^2 - 1)/(A*z-1)^2*dz;
a = [real(w)  real(w+dw)]; b = [imag(w) imag(w+dw)];
line(a,b, 'linewidth', 3, 'color', color2);
end