%   Author: Reija Ristilä <reija@hut.fi>
%   Last modified: 1999-03-25 19:33:06 EET

% Ex 9.2

% coordinates X,Y
r=linspace(0,1,20);
theta=linspace(0,2*pi,100);
X = r'*cos(theta);
Y = r'*sin(theta);

% 1/ prob 2
aa=[0.9 0.5 -0.1 0.05]; 
for i=1:4
  a=aa(i); b=1-a^2;
  num = (1-X.^2-Y.^2);
  denom = 2*pi*((X-a).^2+(Y-b).^2);
  pk = num./denom;
  subplot(2,2,i),surf(X,Y,pk),shading flat, axis([-1 1 -1 1 0 5]);
end

% R & TH
R = r'*ones(size(theta));
TH = ones(size(r'))*theta;

% 2/ prob 3

nn=[1 2 3 4];
for i=1:4
  n=nn(i);
  subplot(2,2,i),mesh(X,Y,R.^n.*cos(n*TH));
end

for i=1:4
  n=nn(i);
  subplot(2,2,i),mesh(X,Y,R.^n.*sin(n*TH));
end

% 3/ prob 4
%z=zeros(19,63);
%for j=1:19
%  R=r(j);
%  for k=1:63
%    TH=t(k);
%    a=(1+R)./(2*R.^2-2*R).*(-2*R.*(R-1)./(R+1)+(R-1).^2.*sin(TH));
%    %b=(1+R)./(2*R.^2-2*R).*(2*(R-1).^2.*cos(TH));
%    z(j,k)=a;
%  end
%end
