Hi Aisha,
Taking the second part first, a phase angle can be either positive or negative. The Matlab angle function uses the usual convention, -pi < theta <= pi. But phase angle has a 2*pi ambiguity, in that you can add 2*pi to any phase angle without changing the essential results. So if you want all positive angles, you can change the range to 0 <= theta < 2*pi, by using the mod function for example:
mod(angle(z),2*pi)
This process leaves the previous positive angles alone and adds 2*pi to the negative ones.
Suppose you have the function u = z in the complex plane. As you go around the origin, starting on the negative x axis, Its phase angle varies from -pi to pi. Using the parula colorbar, the plot on the upper left shows this, with a 2*pi discontinuity from yellow to purple on the negative x axis. That discontinuity is not really consequential, though, so it helps to use a colorbar that goes through the colors and ends where it began. Then the 2*pi disontinuity is not visible. The upper right plot below does that, using a colormap called wheelmap (see details there). Wheelmap is best used in conjuction with the command
caxis([-pi pi])
so the colorbar always goes between those limits. (Otherwise autoscaling of the limits can mess things up).
The function u = sqrt(z) is also shown below. In that case, though, the discontinuity on the negative x axis is real, since the size of the discontinuity is pi and not 2*pi.
Magnitude plots are not as interesting as angle plots. The below show magnitude (actually log10 magnitude) and phase for a quadratic function and for a sine function. The roots of the quadratic are quite visible. You can do the same thing as with function 1 in the code by just plugging in your function
u = (0.5) +(0.5+0.8660i)*z +(-0.25+0.4330i)*z.^2
It turns out that the zeros of the function are at
-0.4967 + 0.8718i and -0.5033 + 0.8603i
though, so you will have to change xx and yy in the code below to create much smaller upper and lower limits and much finer spacing, in order to resolve the zeros.
% create complex plane
xx = -6:.01:6;
yy = -6:.01:6;
[x y] =meshgrid(xx,yy);
z = x+i*y;
% z and sqrt(z)
figure(1)
u = z;
subplot(2,2,1)
surf(x,y,angle(u),'edgecolor','none')
view(2)
colormap(gca,'parula');
colorbar
title('u = z, parula')
subplot(2,2,2)
surf(x,y,angle(u),'edgecolor','none')
caxis([-pi pi])
view(2)
colormap(gca,wheelmap)
colorbar
title('u = z, wheelmap')
u = sqrt(z);
subplot(2,2,3)
surf(x,y,angle(u),'edgecolor','none')
view(2)
colormap(gca,wheelmap(10))
caxis([-pi pi])
colorbar
title('u = sqrt(z), wheelmap')
% polynomial and sine
figure(2)
% function 1
u = (z-(2+2i)).*(z-(-2+i));
subplot(2,2,1)
surf(x,y,log10(abs(u)),'edgecolor','none')
view(2)
colormap(gca,parula)
colorbar
title('polynomial, mag')
subplot(2,2,2)
surf(x,y,angle(u),'edgecolor','none')
view(2)
colormap(gca,wheelmap)
caxis([-pi pi])
colorbar
title('polynomial, angle')
% function 2
u = sin(z);
subplot(2,2,3)
surf(x,y,log10(abs(u)),'edgecolor','none')
view(2)
colormap(gca,parula)
colorbar
title('sin(z), mag')
subplot(2,2,4)
surf(x,y,angle(u),'edgecolor','none')
view(2)
colormap(gca,wheelmap)
caxis([-pi pi])
colorbar
title('sin(z), angle')
.
function map = wheelmap(n,sat,rev)
% 100x3 colormap that goes in color wheel order. Uses hsv2rgb.
% colors can be circularly shifted, and order can be reversed, see below.
% bottom (row 1) to top (row 100) of colorbar is
% red -> yellow -> green -> cyan -> blue -> magenta->red
%
% optional inputs:
% integer n produces a circular shift in the colors.
% Positive vals circularly shift colors up on the colorbar, negative vals
% shift them down. Default is 0.
% Values of n for particular colors at bottom and top:
% n = 0 red->red
% -17 yellow->yellow
% -33 green->green
% 50 cyan->cyan
% 33 blue->blue
% 17 magenta->magenta
%
% sat is the hsv saturation. Default is 2/3.
%
% if rev is 'r' (in quotes) then colors go in reverse order.
% To bypass unneeded options use [], for example wheelmap([],[],'r')
%
% for angle(z) in the complex plane, the 2pi discontinuity does not appear
% and the colors go around z = 0 in counterclockwise order.
% for a general complex function, use
% caxis([-pi pi])
% to bypass autoscaling. Then branch cut discontinuities will be correct.
if nargin == 0
n = 0; sat = 2/3; rev = 'x';
elseif nargin == 1
sat = 2/3; rev = 'x';
elseif nargin == 2
rev = 'x'
end
if isempty(n)
n = 0;
end
if isempty(sat)
sat = 2/3;
end
h = linspace(0,1,100)';
s = sat*ones(size(h));
v = ones(size(h));
map = hsv2rgb([h s v]);
map = circshift(map,n,1);
if rev == 'r'
map = flipud(map);
end
