% June 2014, P.G. Martinsson
%
% This program compares different techniques for evaluating the rank of
% interaction between "source" and "target" regions of different
% geometries.
%
% This has not been proofed at all - use at your own risk!

function main_helmholtz_expansions

%LOCAL_simple
%LOCAL_box_to_circle
%LOCAL_box_to_line
LOCAL_box_to_box
%LOCAL_box_to_circle_quad
%LOCAL_box_to_line_quad
%LOCAL_box_to_line_randomized

return

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

function LOCAL_simple

%xt = [3,2];
xt = [1,sqrt(3)];
rt = norm(xt);
tt = atan2(xt(2),xt(1));

xs = [1,1];
rs = norm(xs);
ts = atan2(xs(2),xs(1));

P = 25;
f_approx = 0;
for p = -P:P
  f_approx = f_approx + besselj(p,rs)*besselh(p,rt)*exp(-1i*p*(tt - ts));
end

f = besselh(0,norm(xt - xs));
fprintf(1,'Approximation = %15.12f + i*%15.12f\n',real(f_approx),imag(f_approx))
fprintf(1,'Exact         = %15.12f + i*%15.12f\n',real(f),imag(f))

return

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

function LOCAL_box_to_circle

rng(0)

NS     = 300;       % Number of source points.
NT     = 200;       % Number of target points.
kh     = 5.5;       % Wave number.
lambda = 2*pi/kh;   % Wave length.
cc     = [0;0];     % Center of expansion.
k_max  = 80;        % Max rank to be investigated.

%%% Set up the source geometry
xxs = 0.5 - rand(2,NS);

%%% Set up the target geometry
tt  = linspace(0,2*pi*(1 - 1/NT),NT);
xxt = 1.5*[cos(tt);sin(tt)];

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(0.5*[-1,1,1,-1,-1],0.5*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',20)
set(han(3),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry: Source region has side = %0.3f lambda',1/lambda))

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
RR      = sqrt(DD1.*DD1 + DD2.*DD2);
A       = besselh(0,kh*RR);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(xxs,kh,cc,P_max);
BB        = LOCAL_T_tfo(xxt,kh,cc,P_max);
for p = 0:P_max
  ind                      = (P_max+1-p):(P_max+1+p);
  err_mpole([2*p+2,2*p+3]) = norm(A - BB(:,ind)*CC(ind,:));
end

figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors: Source region has side = %0.3f lambda',1/lambda))
xlabel('k')
ylabel('||A - A_k||/||A||')
keyboard
    
return

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

function LOCAL_box_to_line

rng(0)

NS     = 300;       % Number of source points.
NT     = 200;       % Number of target points.
kh     = 5.5;       % Wave number.
lambda = 2*pi/kh;   % Wave length.
cc     = [0;0];     % Center of expansion.
k_max  = 80;        % Max rank to be investigated.

%%% Set up the source geometry
xxs = 0.5 - rand(2,NS);

%%% Set up the target geometry
nside = ceil(NT/4);
h     = 3/nside;
tt    = linspace(-1.5,1.5-h,nside);
xxt   = [tt,1.5*ones(1,nside),-tt,-1.5*ones(1,nside);...
         -1.5*ones(1,nside),tt,1.5*ones(1,nside),-tt];

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(0.5*[-1,1,1,-1,-1],0.5*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',20)
set(han(3),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry: Source region has side = %0.3f lambda',1/lambda))

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
RR      = sqrt(DD1.*DD1 + DD2.*DD2);
A       = besselh(0,kh*RR);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(xxs,kh,cc,P_max);
BB        = LOCAL_T_tfo(xxt,kh,cc,P_max);
for p = 0:P_max
  ind                      = (P_max+1-p):(P_max+1+p);
  err_mpole([2*p+2,2*p+3]) = norm(A - BB(:,ind)*CC(ind,:));
end

figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors: Source region has side = %0.3f lambda',1/lambda))
xlabel('k')
ylabel('||A - A_k||/||A||')
keyboard
    
return

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

function LOCAL_box_to_box

rng(0)

NS     = 300;       % Number of source points.
NT     = 300;       % Number of target points.
kh     = 5.5;       % Wave number.
lambda = 2*pi/kh;   % Wave length.
cc     = [0;0];     % Center of expansion.
k_max  = 80;        % Max rank to be investigated.

%%% Set up the source geometry
xxs = 0.5 - rand(2,NS);

%%% Set up the target geometry
xxt = [2.5 - rand(1,NT);0.5 - rand(1,NT)];

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(0.5*[-1,1,1,-1,-1],0.5*[-1,-1,1,1,-1],'k',...
           2 + 0.5*[-1,1,1,-1,-1],0.5*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',20)
set(han(4),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry: Source region has side = %0.3f lambda',1/lambda))

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
RR      = sqrt(DD1.*DD1 + DD2.*DD2);
A       = besselh(0,kh*RR);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(xxs,kh,cc,P_max);
BB        = LOCAL_T_tfo(xxt,kh,cc,P_max);
for p = 0:P_max
  ind                      = (P_max+1-p):(P_max+1+p);
  err_mpole([2*p+2,2*p+3]) = norm(A - BB(:,ind)*CC(ind,:));
end

figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors: Source region has side = %0.3f lambda',1/lambda))
xlabel('k')
ylabel('||A - A_k||/||A||')
keyboard
    
return

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

function LOCAL_box_to_circle_quad

rng(0)

kh     = 5.5;   % Helmholtz parameter
cc     = [0;0]; % Center of geometry.
rad    = 1.5;   % Target circle has radius "rad".
a_in   = 0.5;   % Inner box is of size 2*a_in x 2*a_in.
NT     = 200;   % Number of points along the outer contour.
n      = 20;    % Number of Chebyshev nodes along a side of inner box.
lambda = 2*pi/kh;   % Wave length.
k_max  = 80;        % Max rank to be investigated.

%%% Set up the source geometry.
[x,w]      = LOCAL_clencurt(n,-a_in,a_in);
x          = x(end:(-1):1)';
w          = w';
[XX1,XX2]  = meshgrid(x);
xxs        = [reshape(XX1,1,n*n);...
              reshape(XX2,1,n*n)];
wws        = reshape(w' * w,1,n*n);
NS         = size(xxs,2);
fprintf(1,'Source box has side        = %7.3f lambda     ppw = %0.3f\n',...
        2*a_in/lambda,n*lambda/(2*a_in))

%%% Set up the target geometry.
tt   = linspace(0,2*pi*(1 - 1/NT),NT);
xxt  = rad*[cos(tt);sin(tt)];
wwt  = (rad*2*pi/NT)*ones(1,NT);
circ = 2*pi*rad/lambda; % Circumference measured in wavelengths.
fprintf(1,'Target circle has diameter = %7.3f lambda     ppw = %0.3f\n',...
        2*rad/lambda,NT/circ)

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
RR      = sqrt(DD1.*DD1 + DD2.*DD2);
A       = besselh(0,kh*RR).*sqrt(wwt' * wws);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(xxs,kh,cc,P_max);
BB        = LOCAL_T_tfo(xxt,kh,cc,P_max);
for p = 0:P_max
  J                        = (P_max+1-p):(P_max+1+p);
  A_approx                 = (BB(:,J)*CC(J,:)).*sqrt(wwt' * wws);
  err_mpole([2*p+2,2*p+3]) = norm(A - A_approx);
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',20)
set(han(3),'MarkerSize',20)
axis equal
axis off
title(sprintf('The geometry: Source region has side = %0.3f lambda',2*a_in/lambda))

%%% Plot the errors.
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors: Source region has side = %0.3f lambda',1/lambda))
xlabel('k')
ylabel('||A - A_k||/||A||')

%%% Plot a skeleton
figure(3); subplot(1,1,1); hold off;
acc  = 1e-10;
k_id = max(find(err_ID > (ss(1)*acc)));
han = plot(xxs(1,:),xxs(2,:),'k.',...
           a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,ind(1:k_id)),xxs(2,ind(1:k_id)),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'MarkerSize',20)
set(han(1),'MarkerEdgeColor',0.7*[1,1,1])
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',20)
set(han(4),'MarkerSize',20)
axis equal
axis off
title(sprintf('Skeleton points: eps=%0.1e  k=%d   side of source box = %0.3f lambda',acc,k_id,2*a_in/lambda))
keyboard
return

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

function LOCAL_box_to_line_quad

rng(0)

kh     = 51;   % Helmholtz parameter
cc     = [0;0]; % Center of geometry.
a_in   = 0.5;   % Inner box is of size 2*a_in x 2*a_in.
a_out  = 1.5;   % Outer box is of size 2*a_out x 2*a_out.
n      = 30;    % Number of Chebyshev nodes along a side of inner box.
n_out  = 200;    % Number of Chebyshev nodes along a side of outer box.
lambda = 2*pi/kh;   % Wave length.
k_max  = 150;        % Max rank to be investigated.

%%% Set up the source geometry.
[x,w]      = LOCAL_clencurt(n,-a_in,a_in);
x          = x(end:(-1):1)';
w          = w';
[XX1,XX2]  = meshgrid(x);
xxs        = [reshape(XX1,1,n*n);...
              reshape(XX2,1,n*n)];
wws        = reshape(w' * w,1,n*n);
NS         = size(xxs,2);
fprintf(1,'Source box has side = %7.3f lambda     ppw = %0.3f\n',...
        2*a_in/lambda,n*lambda/(2*a_in))

%%% Set up the target geometry.
[x,w] = LOCAL_clencurt(n_out,-a_out,a_out);
x     = x(end:(-1):1)';
w     = w';
xxt   = [x,a_out*ones(1,n_out-2),x(n_out:(-1):1),-a_out*ones(1,n_out-2);...
         -a_out*ones(1,n_out-1),x,a_out*ones(1,n_out-2),x(n_out:(-1):2)];
wwt   = [2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1))];
NT    = size(xxt,2);
fprintf(1,'Source box has side = %7.3f lambda     ppw = %0.3f\n',...
        2*a_out/lambda,n_out*lambda/(2*a_out))

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
RR      = sqrt(DD1.*DD1 + DD2.*DD2);
A       = besselh(0,kh*RR).*sqrt(wwt' * wws);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the ID
[Q,R,ind] = LOCAL_pivoted_QR3(A,k_max);
err_ID    = [ss(1),NaN*ones(1,k_max)];
for i = 2:length(kk)
  k         = kk(i);
  err_ID(i) = norm(A(:,ind) - Q(:,1:k)*R(1:k,:));
end

%%% Perform the multipole expansions.
err_mpole = [ss(1),NaN*ones(1,k_max)];
P_max     = floor(k_max/2)-1;
CC        = LOCAL_T_ofs(xxs,kh,cc,P_max);
BB        = LOCAL_T_tfo(xxt,kh,cc,P_max);
for p = 0:P_max
  J                        = (P_max+1-p):(P_max+1+p);
  A_approx                 = (BB(:,J)*CC(J,:)).*sqrt(wwt' * wws);
  err_mpole([2*p+2,2*p+3]) = norm(A - A_approx);
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',15)
set(han(3),'MarkerSize',15)
axis equal
axis off
title(sprintf('The geometry: Source region has side = %0.3f lambda',2*a_in/lambda))

%%% Plot a skeleton
figure(1); subplot(1,1,1); hold off;
acc  = 1e-12;
k_id = max(find(err_ID > (ss(1)*acc)));
han = plot(xxs(1,:),xxs(2,:),'k.',...
           a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,ind(1:k_id)),xxs(2,ind(1:k_id)),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'MarkerSize',15)
set(han(1),'MarkerEdgeColor',0.7*[1,1,1])
set(han(2),'LineWidth',2)
set(han(3),'MarkerSize',15)
set(han(4),'MarkerSize',15)
axis equal
axis off
title(sprintf('Skeleton points: eps=%0.1e  k=%d   side of source box = %0.3f lambda',acc,k_id,2*a_in/lambda))

%%% Plot the errors.
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'x-',...
         kk,(1/ss(1))*err_ID,'x-',...
         kk,(1/ss(1))*err_mpole,'x-');
legend('svd','ID','mpole')
title(sprintf('Expansion errors: Source region has side = %0.3f lambda',1/lambda))
xlabel('k')
ylabel('||A - A_k||/||A||')

keyboard
return

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

function LOCAL_box_to_line_randomized

kh     = 51;   % Helmholtz parameter
%kh     = 5;   % Helmholtz parameter
cc     = [0;0]; % Center of geometry.
a_in   = 0.5;   % Inner box is of size 2*a_in x 2*a_in.
a_out  = 1.5;   % Outer box is of size 2*a_out x 2*a_out.
n      = 30;    % Number of Chebyshev nodes along a side of inner box.
n_out  = 100;    % Number of Chebyshev nodes along a side of outer box.
lambda = 2*pi/kh;   % Wave length.
k_max  = 150;        % Max rank to be investigated.

%%% Set up the source geometry.
[x,w]      = LOCAL_clencurt(n,-a_in,a_in);
x          = x(end:(-1):1)';
w          = w';
[XX1,XX2]  = meshgrid(x);
xxs        = [reshape(XX1,1,n*n);...
              reshape(XX2,1,n*n)];
wws        = reshape(w' * w,1,n*n);
NS         = size(xxs,2);
fprintf(1,'Source box has side = %7.3f lambda     ppw = %0.3f\n',...
        2*a_in/lambda,n*lambda/(2*a_in))

%%% Set up the target geometry.
[x,w] = LOCAL_clencurt(n_out,-a_out,a_out);
x     = x(end:(-1):1)';
w     = w';
xxt   = [x,a_out*ones(1,n_out-2),x(n_out:(-1):1),-a_out*ones(1,n_out-2);...
         -a_out*ones(1,n_out-1),x,a_out*ones(1,n_out-2),x(n_out:(-1):2)];
wwt   = [2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1)),2*w(1),w(2:(n_out-1))];
NT    = size(xxt,2);
fprintf(1,'Source box has side = %7.3f lambda     ppw = %0.3f\n',...
        2*a_out/lambda,n_out*lambda/(2*a_out))

%%% Construct the matrix A.
DD1     = xxt(1,:)'*ones(1,NS) - ones(NT,1)*xxs(1,:);
DD2     = xxt(2,:)'*ones(1,NS) - ones(NT,1)*xxs(2,:);
RR      = sqrt(DD1.*DD1 + DD2.*DD2);
A       = besselh(0,kh*RR).*sqrt(wwt' * wws);
ss      = svd(A);
kk      = 0:k_max;
err_svd = ss(1:(k_max+1));
fprintf(1,'||A|| = %0.10f\n',ss(1))

%%% Perform the randomized sampling.
n_exp = 1; % Number of experiments.
err_rand = [ss(1),NaN*ones(1,k_max)];
ERR = err_rand'*ones(1,n_exp);
for j = 1:size(ERR,2)
  disp(j)
  Omega = randn(NS,k_max);
  Y     = A*Omega;
  for i = 2:length(kk)
    k = kk(i);
    [Q,~] = qr(Y(:,1:k),0);
    ERR(i,j) = norm(A - Q*(Q'*A));
  end
end

%%% Plot the geometry.
figure(1); subplot(1,1,1); hold off;
han = plot(a_in*[-1,1,1,-1,-1],a_in*[-1,-1,1,1,-1],'k',...
           xxs(1,:),xxs(2,:),'r.',...
           xxt(1,:),xxt(2,:),'b.');
set(han(1),'LineWidth',2)
set(han(2),'MarkerSize',15)
set(han(3),'MarkerSize',15)
axis equal
axis off
title(sprintf('The geometry: Source region has side = %0.3f lambda',2*a_in/lambda))

%%% Plot the errors.
figure(2); subplot(1,1,1); hold off;
semilogy(kk,(1/ss(1))*err_svd,'bx-',...
         kk,(1/ss(1))*ERR,'rx-');
legend('svd','randomized')
%title(sprintf('Expansion errors: Source region has side = %0.3f lambda',1/lambda))
xlabel('k')
ylabel('||A - A_k||/||A||')

keyboard
return

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

function T = LOCAL_T_ofs(xxs,kh,cc,P)

dd1 = xxs(1,:) - cc(1);
dd2 = xxs(2,:) - cc(2);
rr  = sqrt(dd1.*dd1 + dd2.*dd2);
tt  = atan2(dd2,dd1);
NS  = size(xxs,2);
T   = zeros(2*P+1,NS);

for p = (-P):P
  T(p+P+1,:) = besselj(p,kh*rr).*exp(1i*p*tt);
end

return

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

function T = LOCAL_T_tfo(xxt,kh,cc,P)

dd1 = (xxt(1,:) - cc(1))';
dd2 = (xxt(2,:) - cc(2))';
rr  = sqrt(dd1.*dd1 + dd2.*dd2);
tt  = atan2(dd2,dd1);
NT  = size(xxt,2);
T   = zeros(NT,2*P+1);

for p = (-P):P
  T(:,p+P+1) = besselh(p,kh*rr).*exp(-1i*p*tt);
end

return

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

function [Q,R,ind] = LOCAL_pivoted_QR3(A,k)

% This function orthogonalizes the COLUMNS of A
% It uses a modified Gram-Schmidt with column pivoting

m = size(A,1);
n = size(A,2);

R = zeros(min(m,n),n);
Q = A;
ind = 1:n;

for j = 1:k
    [trash, j_max] = max(sum(Q(:,j:n).*Q(:,j:n)));
    j_max          = j_max + j - 1;
    Q(:,[j,j_max]) = Q(:,[j_max,j]);
    R(:,[j,j_max]) = R(:,[j_max,j]);
    ind([j,j_max]) = ind([j_max,j]);
    r_jj   = norm(Q(:,j));
    Q(:,j) = Q(:,j)/r_jj;
    Q(:,j) = Q(:,j) - Q(:,1:(j-1))*(Q(:,1:(j-1))'*Q(:,j));
    Q(:,j) = Q(:,j)/norm(Q(:,j));
%    Q(:,j) = Q(:,j) - Q(:,1:(j-1))*(Q(:,1:(j-1))'*Q(:,j));
%    Q(:,j) = Q(:,j)/norm(Q(:,j));
    R(j,j) = r_jj;
    rr     = Q(:,j)'*Q(:,(j+1):n);
    R(j,(j+1):n) = rr;
    Q(:,(j+1):n) = Q(:,(j+1):n) - Q(:,j)*rr;
end

Q = Q(:, 1:min(m,n));

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the Chebyshev nodes on [-1,1].
% It also computes a differentiation operator.
% It is modified from a code by Nick Trefethen.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [D,x] = LOCAL_cheb(N)

if N==0
  D=0; 
  x=1; 
  return
end
x = cos(pi*(0:N)/N)'; 
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';                  
D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
D  = D - diag(sum(D'));                 % diagonal entries

return

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

function [x,w] = LOCAL_clencurt(N1,a,b)

% Adapted from "fclencurt.m" by Greg von Winckel.
N=N1-1; 
bma=b-a;
c=zeros(N1,2);
c(1:2:N1,1)=(2./[1 1-(2:2:N).^2 ])'; 
c(2,2)=1;
f=real(ifft([c(1:N1,:);c(N:-1:2,:)]));
w=bma*([f(1,1); 2*f(2:N,1); f(N1,1)])/2;
x=0.5*((b+a)+N*bma*f(1:N1,2));

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The outgoing-from-sources translation operator:
%   T(  1,j) = 1;
%   T(p+1,j) = -(1/p)*(xxj - cc)^p
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

