%                          Program dirch
%
%  This program solves the Dirichlet problem in the unit disk for given
%  boundary data f. The program calculates the approximate Fourier
%  coefficients of g using the Matlab fast fourier transform and then
%  sums the terms n = 0,. . . , 63. User must provide an array smart 
%  function mfile called f.m for the boundary data.  f must be periodic
%  of period 2pi.



N= 128;
M = 21;

r = 0:.05:1;
theta= 0: 2*pi/N :2*pi;
X = r'* cos(theta);
Y = r'* sin(theta);


% Calculate the Fourier coefficients of the boundary data f 


ff = f(theta(1:N));

dd = fft(ff)/N;

d0 = dd(1);
d = dd(2:N);


U = d0*ones(size(X));

for k = 1:N/2 -1
     U = U +(r.^k)'*(d(k)*exp(i*k*theta) + d(N-k)*exp(-i*k*theta));
end
     
surf(X,Y,U); shading flat

m = min(min(U))-.5;
hold on

plot3([0 1], [0 0], [m m], 'r')
plot3([.9 1], [-.1 0], [m m], 'r')
plot3([.9 1], [ .1 0], [m m], 'r')
plot3([0 0], [0 1], [m m], 'r')
plot3([-.1 0], [.9 1], [m m], 'r')
plot3([ .1 0], [.9 1], [m m], 'r')
plot3( cos(theta), sin(theta), m*ones(size(theta)), 'r' )


view(135,30)
hold off