function main04

%LOCAL_single_bump

%LOCAL_two_bumps

%LOCAL_src_and_trg

LOCAL_src_and_trg_helm

%LOCAL_circle

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_single_bump

nside = 100;
a     = 3;
h     = 2*a/(nside-1);
xxc   = [0,0];
r     = 1;
p     = 2;
[XX1,XX2] = meshgrid(linspace(-a,a,nside));
xx1   = reshape(XX1,nside*nside,1);
xx2   = reshape(XX2,nside*nside,1);

%%% Construct the analytic solution.
[u,dudx1,dudx2,lapu] = LOCAL_soln(xx1,xx2,xxc,r);

%%% Do numerical differentiation.
cc = LOCAL_1D_stencil(p);
STENCIL = [[0;0;2*cc(1)],zeros(3,4*p)];
for j = 1:p
  STENCIL(:,1 + 4*(j-1) + (1:4)) = [ 0,            0,      j,     -j;...
                                     j,           -j,      0,      0;...
                                     cc(j+1),cc(j+1),cc(j+1),cc(j+1)];
end
A = (1/(h*h))*LOCAL_regular_stencil(nside,nside,STENCIL);
g = A*u;
I_ext = find( ((a - abs(xx1)) < (p-0.5)*h) | ...
             ((a - abs(xx2)) < (p-0.5)*h) );
g(I_ext) = NaN*ones(length(I_ext),1);

figure(1)
subplot(1,2,1)
mesh(XX1,XX2,reshape(lapu,nside,nside))
axis equal
title('The body load g.')
subplot(1,2,2)
mesh(XX1,XX2,reshape(u,nside,nside))
axis tight
title('The solution u.')
keyboard
export_fig fig_laplace_soln1.pdf -pdf -transparent

figure(1)
subplot(1,2,1)
mesh(XX1,XX2,reshape(dudx1,nside,nside))
axis([-a,a,-a,a,-1,1])
title('du/dx1')
subplot(1,2,2)
mesh(XX1,XX2,reshape(dudx2,nside,nside))
axis([-a,a,-a,a,-1,1])
title('du/dx2')
export_fig fig_laplace_soln2.pdf -pdf -transparent

figure(3)
subplot(1,1,1)
mesh(XX1,XX2,log10(abs(reshape(g - lapu,nside,nside)) + 1e-15))
title('Discrepancy between exact g and computed g.')
zlabel('10-logarithm of discrepancy')

keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_two_bumps

nside = 100;
a     = 3;
h     = 2*a/(nside-1);
p     = 1;
[XX1,XX2] = meshgrid(linspace(-a,a,nside));
xx1   = reshape(XX1,nside*nside,1);
xx2   = reshape(XX2,nside*nside,1);

%%% Construct the analytic solution.
[v,dvdx1,dvdx2,lapv] = LOCAL_soln(xx1,xx2,[ 1,0],0.5);
[w,dwdx1,dwdx2,lapw] = LOCAL_soln(xx1,xx2,[-1,0],0.5);
u     = v - w;
dudx1 = dvdx1 - dwdx1;
dudx2 = dvdx2 - dwdx2;
lapu  = lapv - lapw;

%%% Do numerical differentiation.
cc = LOCAL_1D_stencil(p);
STENCIL = [[0;0;2*cc(1)],zeros(3,4*p)];
for j = 1:p
  STENCIL(:,1 + 4*(j-1) + (1:4)) = [ 0,            0,      j,     -j;...
                                     j,           -j,      0,      0;...
                                     cc(j+1),cc(j+1),cc(j+1),cc(j+1)];
end
A = (1/(h*h))*LOCAL_regular_stencil(nside,nside,STENCIL);
f = A*u;
I_ext = find( ((a - abs(xx1)) < (p-0.5)*h) | ...
              ((a - abs(xx2)) < (p-0.5)*h) );
f(I_ext) = NaN*ones(length(I_ext),1);

figure(1)
subplot(1,2,1)
mesh(XX1,XX2,reshape(lapu,nside,nside))
title('The right-hand side f')
subplot(1,2,2)
mesh(XX1,XX2,reshape(u,nside,nside))
title('The solution u')

figure(2)
subplot(1,2,1)
mesh(XX1,XX2,reshape(dudx1,nside,nside))
title('du/dx1')
subplot(1,2,2)
mesh(XX1,XX2,reshape(dudx2,nside,nside))
title('du/dx2')

figure(3)
subplot(1,1,1)
mesh(XX1,XX2,log10(abs(reshape(f - lapu,nside,nside)) + 1e-15))
title('Discrepancy between exact f and computed f.')
zlabel('10-logarithm of discrepancy')

keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_src_and_trg

rng(5)
ns    = 5;
nside = 20;

xxs   = rand(2,ns);
%qq    = ones(ns,1);
%qq    = -(-1).^(1:ns);
%qq    = qq' - mean(qq);
trash = randn(150,1);
qq    = randn(ns,1);
qq    = qq - mean(qq);

[XX1,XX2] = meshgrid(linspace(-1,4,5*nside+1),linspace(-1,2,3*nside+1));
nt        = numel(XX1);

xxt = [reshape(XX1,1,nt);...
       reshape(XX2,1,nt)];

dd1 = xxt(1,:)'*ones(1,ns) - ones(nt,1)*xxs(1,:);
dd2 = xxt(2,:)'*ones(1,ns) - ones(nt,1)*xxs(2,:);
uu  = -0.5*log(dd1.*dd1 + dd2.*dd2)*qq;
UU  = reshape(uu,size(XX1,1),size(XX1,2));

%figure(1)
%hold off
%plot([0,1,1,0,0],[0,0,1,1,0],'r',...
%     [2,3,3,2,2],[0,0,1,1,0],'b','LineWidth',3)
%axis equal
%axis([-0.1,3.1,-0.1,1.1])
%axis off
%keyboard
%%export_fig laplace_src_trg_fig0.pdf -pdf -transparent

figure(1)
hold off
han = plot3([0,1,1,0,0],[0,0,1,1,0],0*ones(1,5),'r',...
            [2,3,3,2,2],[0,0,1,1,0],0*ones(1,5),'b',...
            xxs(1,:),xxs(2,:),0*ones(1,ns),'r.');
hold on
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(2),'MarkerSize',20)
axis([-1,4,-1,2,-3,3])
view([20,30])
keyboard
%export_fig laplace_src_trg_fig1.pdf -pdf -transparent

figure(1)
hold off
han = plot3([0,1,1,0,0],[0,0,1,1,0],0*ones(1,5),'r',...
            [2,3,3,2,2],[0,0,1,1,0],0*ones(1,5),'b',...
            xxs(1,:),xxs(2,:),0*ones(1,ns),'r.');
hold on
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(2),'MarkerSize',20)
mesh(XX1,XX2,UU);
axis([-1,4,-1,2,-3,3])
view([20,30])
keyboard
%export_fig laplace_src_trg_fig2.pdf -pdf -transparent

figure(1)
h = 1/nside;
ind = find( (xxt(1,:) > (2-0.5*h)) & ...
            (xxt(1,:) < (3+0.5*h)) & ...
            (xxt(2,:) > (0-0.5*h)) & ...
            (xxt(2,:) < (1+0.5*h)) );
VV = reshape(uu(ind),nside+1,nside+1);
YY1 = reshape(xxt(1,ind),nside+1,nside+1);
YY2 = reshape(xxt(2,ind),nside+1,nside+1);
hold off
han = plot3([0,1,1,0,0],[0,0,1,1,0],0*ones(1,5),'r',...
            [2,3,3,2,2],[0,0,1,1,0],0*ones(1,5),'b',...
            xxs(1,:),xxs(2,:),0*ones(1,ns),'r.');
hold on
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(2),'MarkerSize',20)
han = mesh(XX1,XX2,UU);
set(han,'EdgeColor',0.7*[1,1,1])
mesh(YY1,YY2,VV)
axis([-1,4,-1,2,-3,3])
view([20,30])
keyboard
%export_fig laplace_src_trg_fig3.pdf -pdf -transparent

figure(1)
hold off
han = plot3([0,1,1,0,0],[0,0,1,1,0],0*ones(1,5),'r',...
            [2,3,3,2,2],[0,0,1,1,0],0*ones(1,5),'b',...
            xxs(1,:),xxs(2,:),0*ones(1,ns),'r.');
hold on
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(2),'MarkerSize',20)
mesh(YY1,YY2,VV)
axis([-1,4,-1,2,min([uu(ind);-0.2]),max([uu(ind);0.2])])
view([20,30])
keyboard
export_fig laplace_src_trg_fig4bb.pdf -pdf -transparent

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_src_and_trg_helm

rng(5)
ns     = 5;
kh     = 40;
lambda = 2*pi/kh;
nside  = 30;
fprintf(1,'The boxes has side length %0.3f lambda\n',1/lambda)

xxs   = rand(2,ns);
%qq    = ones(ns,1);
%qq    = -(-1).^(1:ns);
%qq    = qq' - mean(qq);
qq    = randn(ns,1) + 1i*randn(ns,1);
qq    = qq - mean(qq);

[XX1,XX2] = meshgrid(linspace(-1,4,5*nside+1),linspace(-1,2,3*nside+1));
nt        = numel(XX1);

xxt = [reshape(XX1,1,nt);...
       reshape(XX2,1,nt)];

dd1 = xxt(1,:)'*ones(1,ns) - ones(nt,1)*xxs(1,:);
dd2 = xxt(2,:)'*ones(1,ns) - ones(nt,1)*xxs(2,:);
uu  = besselh(0,kh*sqrt(dd1.*dd1 + dd2.*dd2))*qq;
UU  = reshape(uu,size(XX1,1),size(XX1,2));

%figure(1)
%hold off
%plot([0,1,1,0,0],[0,0,1,1,0],'r',...
%     [2,3,3,2,2],[0,0,1,1,0],'b','LineWidth',3)
%axis equal
%axis([-0.1,3.1,-0.1,1.1])
%axis off
%keyboard
%%export_fig helm_src_trg_fig0.pdf -pdf -transparent

figure(1)
hold off
han = plot3([0,1,1,0,0],[0,0,1,1,0],0*ones(1,5),'r',...
            [2,3,3,2,2],[0,0,1,1,0],0*ones(1,5),'b',...
            xxs(1,:),xxs(2,:),0*ones(1,ns),'r.');
hold on
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(2),'MarkerSize',20)
axis([-1,4,-1,2,-3,3])
title(sprintf('Helmholtz problem. Side of boxes = %0.2f lambda',1/lambda))
view([20,30])
keyboard
export_fig helm_src_trg_fig1high.pdf -pdf -transparent

figure(1)
hold off
han = plot3([0,1,1,0,0],[0,0,1,1,0],0*ones(1,5),'r',...
            [2,3,3,2,2],[0,0,1,1,0],0*ones(1,5),'b',...
            xxs(1,:),xxs(2,:),0*ones(1,ns),'r.');
hold on
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(2),'MarkerSize',20)
mesh(XX1,XX2,real(UU));
axis([-1,4,-1,2,-3,3])
title(sprintf('Helmholtz problem. Side of boxes = %0.2f lambda',1/lambda))
view([20,30])
keyboard
export_fig helm_src_trg_fig2high.pdf -pdf -transparent

figure(1)
h = 1/nside;
ind = find( (xxt(1,:) > (2-0.5*h)) & ...
            (xxt(1,:) < (3+0.5*h)) & ...
            (xxt(2,:) > (0-0.5*h)) & ...
            (xxt(2,:) < (1+0.5*h)) );
VV = reshape(uu(ind),nside+1,nside+1);
YY1 = reshape(xxt(1,ind),nside+1,nside+1);
YY2 = reshape(xxt(2,ind),nside+1,nside+1);
hold off
han = plot3([0,1,1,0,0],[0,0,1,1,0],0*ones(1,5),'r',...
            [2,3,3,2,2],[0,0,1,1,0],0*ones(1,5),'b',...
            xxs(1,:),xxs(2,:),0*ones(1,ns),'r.');
hold on
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(2),'MarkerSize',20)
han = mesh(XX1,XX2,real(UU));
set(han,'EdgeColor',0.7*[1,1,1])
mesh(YY1,YY2,real(VV))
axis([-1,4,-1,2,-3,3])
title(sprintf('Helmholtz problem. Side of boxes = %0.2f lambda',1/lambda))
view([20,30])
keyboard
export_fig helm_src_trg_fig3high.pdf -pdf -transparent

figure(1)
hold off
han = plot3([0,1,1,0,0],[0,0,1,1,0],0*ones(1,5),'r',...
            [2,3,3,2,2],[0,0,1,1,0],0*ones(1,5),'b',...
            xxs(1,:),xxs(2,:),0*ones(1,ns),'r.');
hold on
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(2),'MarkerSize',20)
mesh(YY1,YY2,real(VV))
axis([-1,4,-1,2,min([real(uu(ind));-0.2]),max([real(uu(ind));0.2])])
title(sprintf('Helmholtz problem. Side of boxes = %0.2f lambda',1/lambda))
view([20,30])
keyboard
export_fig helm_src_trg_fig4high.pdf -pdf -transparent

hold on
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(2),'MarkerSize',20)
han = mesh(XX1,XX2,abs(UU));
set(han,'EdgeColor',0.7*[1,1,1])
mesh(YY1,YY2,abs(VV))
axis([-1,4,-1,2,-3,3])
title(sprintf('Helmholtz problem. Side of boxes = %0.2f lambda',1/lambda))
view([20,30])
keyboard
export_fig helm_src_trg_fig5high.pdf -pdf -transparent

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_circle

nt = 100;
nr = 50;
[RR,TT] = meshgrid(linspace(0,0.98,nr),linspace(0,2*pi,nt));
XX1 = RR.*cos(TT);
XX2 = RR.*sin(TT);

P = 200;
U = (1/(2*pi))*ones(size(XX1));
for i = 1:P
  U = U + (1/pi)*(RR.^i).*cos(i*TT);
end

uu = reshape(U,1,nt*nr);
ind = find(uu > 3);
uu(ind) = NaN*ones(size(uu(ind)));
U_trunc = reshape(uu,nt,nr);

figure(1)
hold off
ss = linspace(0,2*pi,500);
plot3(cos(ss),sin(ss),zeros(size(ss)),'k')
hold on
mesh(XX1,XX2,U)
axis equal
axis([-1.1,1.1,-1.1,1.1,-0.1,3])
keyboard
export_fig fig_greenfunction_circle.pdf -pdf -transparent
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [u,dudx1,dudx2,lapu] = LOCAL_soln(XX1,XX2,xxc,r)

YY1     = XX1 - xxc(1);
YY2     = XX2 - xxc(2);
rrsq    = YY1.*YY1 + YY2.*YY2;
mask    = (rrsq < (r*r));
u       = 0.25*mask.*(r*r - rrsq) - 0.25*r*r*(1-mask).*(log(rrsq) - log(r*r));
dudx1   = - 0.5*mask.*YY1 - 0.5*r*r*(1-mask).*YY1./rrsq;
dudx2   = - 0.5*mask.*YY2 - 0.5*r*r*(1-mask).*YY2./rrsq;
lapu    = 1*mask;

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function creates a (2p+1)-stencil for approximating the second 
% derivative of a function. In other words, it creates a vector cc such that
% -u''(0) \approx cc(1) + cc(2)*(u(h)+u(-h)) + cc(3)*(u(2*h)+u(-2*h)) + ...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function cc = LOCAL_1D_stencil(p)

X1 = ones(p,1)*(1:p);
X2 = (2:2:(2*p))' * ones(1,p);
X  = X1.^X2;
cc = zeros(p+1,1);
cc(2:end) = X\[-1;zeros(p-1,1)];
cc(1) = -2*sum(cc(2:end));

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_regular_stencil(nside1,nside2,STENCIL)

nrow    = size(STENCIL,2);
mpad    = max(max(abs(STENCIL(1:2,:))));
mside1  = nside1 + 2*mpad;
mside2  = nside2 + 2*mpad;
iii     = mside2*STENCIL(1,:) + STENCIL(2,:);
jjj     = zeros(1,nrow);
aaa     = STENCIL(3,:);

iiilong = ones(1,nside1*nside2*nrow);
jjjlong = ones(1,nside1*nside2*nrow);
aaalong = ones(1,nside1*nside2*nrow);

for icol = 1:nside1
   ind          = nside2*nrow*(icol-1) + (1:(nside2*nrow));
   offsets      = (mside2+1)*mpad + mside2*(icol-1) + (1:nside2);
   iiilong(ind) = kron(ones(1,length(offsets)),iii) + ...
                  kron(offsets,ones(1,nrow));
   jjjlong(ind) = kron(ones(1,length(offsets)),jjj) + ...
                  kron(offsets,ones(1,nrow));
   aaalong(ind) = kron(ones(1,length(offsets)),aaa);
end
A = sparse(iiilong,jjjlong,aaalong,mside1*mside2,mside1*mside2);

indcore = zeros(1,nside1*nside2);
for icol = 1:nside1
  ind = (icol-1)*nside2 + (1:nside2);
  indcore(ind) = (mpad + icol - 1)*(nside2 + 2*mpad) + mpad + (1:nside2);
end

A = A(indcore,indcore);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

