%
%  w10green.m -- djm -- 08 nov 2010 
%

N = 101;  L = 1;  M = 400;

x0 = L/3;  y0 = L/4;

ds = L/N;
xx = 0:ds:L;  yy = xx;

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

ug = zeros(size(xg));
%  double sum
for kk = 1:M
    for ll = 1:M
        temp = (sin(kk*pi/L * x0) * sin(ll*pi/L * y0)) / (kk^2 + ll^2);
        ug = ug + temp * sin(kk*pi/L * xg) .* sin(ll*pi/L * yg);
    end
end
   
ug = (-4/pi^2)*ug;

hg = ug - (1/4/pi)*log((xg-x0).^2 + (yg-y0).^2);

%  plots

figure(1);  clf
surf(xx,yy,ug);  shading flat
title('\bf greens function in a square','fontsize',14)
xlabel('\bf x-axis','fontsize',12)
ylabel('\bf y-axis','fontsize',12)
zlabel('\bf u-axis','fontsize',12)

caxis([-0.5 0.1])
axis([0 L 0 L -0.5 0.1])


figure(2);  clf
surf(xx,yy,hg);  shading flat
title('\bf harmonic part of GF in a square','fontsize',14)
xlabel('\bf x-axis','fontsize',12)
ylabel('\bf y-axis','fontsize',12)
zlabel('\bf u-axis','fontsize',12)

byhg = [hg(1,:) hg(end,:) hg(:,1)' hg(:,end)'];

maxhg = max(byhg(:));
minhg = min(byhg(:));

caxis([minhg maxhg])
axis([0 L 0 L minhg maxhg])

hold on
%  boundary values
xb = xx;  yb = yy(1)  +0*yy;
hb = - (1/4/pi)*log((xb-x0).^2 + (yb-y0).^2);
plot3(xb,yb,hb,'k','linewidth',2)

xb = xx;  yb = yy(end)+0*yy;
hb = - (1/4/pi)*log((xb-x0).^2 + (yb-y0).^2);
plot3(xb,yb,hb,'k','linewidth',2)

xb = xx(1)  +0*xx;  yb = yy;
hb = - (1/4/pi)*log((xb-x0).^2 + (yb-y0).^2);
plot3(xb,yb,hb,'k','linewidth',2)

xb = xx(end)+0*xx;  yb = yy;
hb = - (1/4/pi)*log((xb-x0).^2 + (yb-y0).^2);
plot3(xb,yb,hb,'k','linewidth',2)