% code by Colin Macdonald
% comments by Ben Ong
% this code shows us where
% the singularity in the integrand is.

dx = .02;
dy = .02;

xx = -2*pi:dx:2*pi;
yy = -3:dy:4;

[xg,yg] = meshgrid(xx,yy);

Z = xg+1i*yg;

k = .5;
n = 2;

%F = sin(n*Z) .* sin(Z) ./ (1-k*cos(Z));
F = exp(1i*n*Z) .* sin(Z) ./ (1-k*cos(Z)).^2;
Fstr = 'exp(1i*n*Z) .* sin(Z) ./ (1-k*cos(Z)).^2';

%because this is such a "harsh singularity,
% (i.e. around the singularity, the function
%is very large, we set limits on it when we're plotting

figure(1); clf;
F2 = real(F);
LLIM = -500;
ULIM = 500;
pcolor(xg,yg, F2.*(F2<=ULIM & F2>=LLIM) + LLIM.*(F2<LLIM) + ULIM.*(F2>ULIM));
shading flat;
hold on;
[C,H] = contour(xg,yg,real(F),[0 0],'k-');
colorbar
title(['real( ' Fstr ' )']);
%colormap(gray);

figure(2); clf;
subplot(2,1,2);
F2 = imag(F);
max(max(F2))
min(min(F2));
LLIM = -200;
ULIM = 200;
pcolor(xg,yg, F2.*(F2<=ULIM & F2>=LLIM) + LLIM.*(F2<LLIM) + ULIM.*(F2>ULIM));
shading flat;
hold on
contour(xg,yg,imag(F),[0 0], 'k-');
contour(xg,yg,imag(F),[.1 .1], 'r-');
contour(xg,yg,imag(F),[.2 .2], 'g-');
colorbar
title(['imag( ' Fstr ' )']);
%colormap(gray);

% x-section
subplot(2,1,1);
plot(xx,F2(end,:),'k-');