function [z,w,w1] = viscap(eps,V,r,al)

% We compute resonances of D_x^2 + V where V is a step potential and 
% the eigenvalues of  D_x^2 + V - i epsilon x^2. The input is given 
% [ epsilon , V , r ] where V = [ V(1), V(2), ... V(N) ] are values of 
% the potential in  r(K) < x < r(K+1) where r = [r(1), r(2), ... , r(N+1)]. 
% In the output [ z , w, v ], z are approximate eigenvalues of D_x^2 + V - i epsilon x^2, 
% w are the resonances -- but z and w are "pruned" to include only the eigenvalues 
% near 0 and to exclude spurious einvalues with arg z = pi/4. 
% The last output consists of all the computed eigenvalues of 
% the discretization of D_x^2 + V - i epsilon x^2:
% most of them are spurious or due to pseudospectral effect. 

  if ( nargin < 1) 

    eps = 0.5;

end

if   (nargin < 2 ) 

  V = [ 2 -4 1 5];

end

if   (nargin < 3 ) 

  r = [ -2, - 1, 0 , 1 , 2];

end

if (nargin < 4 ) 

  al = pi/2;

end

x=[];

VV = [];

r = [ r(1) - 1 , r , r(length(r))+1];

V = [ 0 , V , 0 ];

l = length(r);

  for j=1:length(r)-1;

y = linspace(r(j),r(j+1),200);

x = [ x , y ];

VV = [VV, V(j)+0*y ]; 

end

y = x;

for j=1:151
	f(j,:)=herm(y,j-1);
end

VVV = zeros(151);

for j=1:151
	for k=1:j
		VVV(k,j)=trapz(y,VV.*f(j,:).*f(k,:));
end
end

for j=1:150
for k=j+1:151
	VVV(k,j)=VVV(j,k);
end
end

A = -hD1(151)^2 + eps*exp(-i*al)*hx2(151) + VVV;

w1 = eig(A);

w1 = w1;

z = prune(w1);

w = respot(V,r);

w = prune(w);

plot(w,'ko','MarkerSize',10,'LineWidth',2)

hold on 

plot(z,'*','MarkerSize',10,'LineWidth',2)

hold off

legend('resonances','viscosity approximation','Location','southwest')
title(['viscosity approximation with \epsilon = ' num2str(eps)])

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

function hD1 = hD1(N)

%return the matrix in the basis of the quantum harmonic oscillator
%of differentiation in x; N is the number of basis elements;

A = zeros(N);

for kk=1:N;
for mm=1:N

k=kk-1;

m=mm-1;

A(kk,mm)=0.5*ha(m)/ha(k)*krd(m+1,k) - 0.5*ha(k)/ha(m)*krd(m,k+1);

%The second option is the derivative acting on hermite polynomials;
%The basis now is hermite polynomials on L2 with measure e^{-|x|^2}dx;

%A(kk,mm)=2*k*ha(k)/ha(m)*krd(k,m+1);

end
end

hD1=A;


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

function hx2 = hx2(N)

%return the matrix in the basis of the quantum harmonic oscillator
%of multiplication by x^2; N is the number of basis elements;

A = zeros(N);

for kk=1:N;
for mm=1:N

	 k=kk-1;
m=mm-1;
A(kk,mm)=0.25*(ha(m)/ha(k)*krd(m+2,k)+ha(k)/ha(m)*krd(m,k+2))+...
  (m+0.5)*krd(m,k);
end
end

hx2=A;

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

function ha = ha(n);

ha = (2^n*factorial(n)*sqrt(pi))^(-1/2);

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

function y = herm(x, n);

y = hermite(x,n);
y = y';

a = (2^n*factorial(n)*sqrt(pi))^(-1/2);

y = a*y.*exp(-x.^2/2);

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

function y = hermite(x, n);

if size(x, 1) == 1;
  x = x'; rot = 1;
end;

ytem = zeros(length(x), 1);
ytem(:,1) = ones(length(x), 1);
ytem(:,2) = 2*x;

if n == 0; y = ytem(:,1); return; end;
if n == 1; y = ytem(:,2); return; end;

if n > 1;

  m = 1;
  p = 2;

  for i = 2:n;
    ynew = 2*x.*ytem(:,p)-2*(i-1)*ytem(:,m);
    tem = m;
    m = p;
    p = tem; 
    ytem(:,p) = ynew;
  end

  y = ynew;

end;

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

function krd = krd(n,m)

  if (n==m) 

    krd=1;

  else

    krd=0;

end

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

function prune = prune(z,R,th)

  if (nargin < 3) 
    th = 0.8;
end

if (nargin < 3)
  R = 120;
end

p =1;

w = [ 1 2 ];

for j=1:length(z)

	if ( real(z(j)) < R ) 

	  w(p)=z(j);

p=p+1;

end
end

z = w;

w = [1,2];

p = 1;

for j=1:length(z)

	rrr =  abs(imag(z(j)))/abs(real(z(j)));

	if ( rrr < th ) 

w(p)=z(j);

p = p+1;

end
end

z = w;

w = [1,2];

p = 1;

for j=1:length(z)

	if (imag(z(j)) > -20) 

	  w(p)=z(j);

p = p+1;

end
end



w = real(w) - i*abs(imag(w));

prune = w;

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

function respot = respot(V,r)

figure

  z = squarepot(V,r);

p=1;

w=[1];

for j=1:length(z)

	  if ( imag(z(j))/(abs(real(z(j))) + 10^-6) < -10 ) 

	    p = p;

	  else 

	    w(p)=z(j);  

p = p+1;

end
end

w = w.^2;

  respot=prune(w);

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

% function l = squarepot(VV, xx, neigs)
%
% Compute resonances of a one-dimensional piecewise-constant potential.
% The value of the potential will be VV(i) on the interval ( xx(i), xx(i+1) ).
%
% Input:
%   VV - Value of the potential on each interval
%   xx - Coordinates of the interval endpoints
%   neigs - Max eigenvalues or resonances to be computed (default: 40)
%
% Output:
%   l - Vector of resolved resonance poles in the lambda (wave number) plane

function l = squarepot(VV,xx, neigs);

if nargin ==0
  VV = [0, 1, 0];
  xx = [-2,-1,1,2];
end

if nargin < 3
  neigs = 40;
end

elt = square_well(xx,VV);
l = checked_resonances(elt,neigs);

clf;
subplot(2,1,1) ;
plot_potential1(VV,xx); 
title('Potential');

subplot(2,1,2);
plot(checked_resonances(elt), '.', 'MarkerSize', 25); 
title('Pole locations');

% elt = square_well(ab, V, h, order)
%
% Discretize a problem where the potential is given by a piecewise constant.
% Input:
%   ab    - breakpoints between regions where the potential is constant
%   V     - potential values on each interval ab(j) to ab(j+1)
%   h     - maximum element size (default: 1)
%   order - maximum element order (default: 20)

function elt = square_well(ab, V, h, order)

if nargin == 0
  ab = [-2, -1, 1, 2];
  V  = [0, -10, 0];
end

if nargin == 1
  V  = [0, ab, 0];
  ab = [-2, -1, 1, 2];
end

if nargin < 3, h = 1;      end
if nargin < 4, order = 20; end

elt = [];
nint = length(V);
i = 1;
for j = 1:nint
  L = ab(j+1)-ab(j);
  nelt = ceil(L/h);
  for k = 1:nelt
    elt(i).a = ab(j) + L*(k-1)/nelt;
    elt(i).b = ab(j) + L*k/nelt;
    elt(i).order = order;
    elt(i).Vtype = 'constant';
    elt(i).V = V(j);
    i = i + 1;
  end
end
% [l] = checked_resonances(elt,neigs,tol)
%
% Compute resonances with two densities in order to check convergence.
% If the same answer occurs to within tol, accept the pole as converged.
% Inputs:
%   elt   - coarse mesh
%   neigs - number of poles desired? (default: 0 -> compute all)
%   tol   - match tolerance (default: 1e-6)

function l = checked_resonances(elt,neigs,tol)

if nargin < 2, neigs = 0;  end
if nargin < 3, tol = 1e-6; end

elt2 = elt;
for j = 1:length(elt)
  elt2(j).order = ceil(elt2(j).order*1.5);
end
l1 = compute_resonances(elt,  neigs);
l2 = compute_resonances(elt2, neigs);

l = compare_eigs(l2,l1,tol);
% [l,V] = compute_resonances(elt, neigs)
%
% Compute resonances via a generalized linear eigenvalue problem.
% Inputs:
%   elt   - Problem description
%   neigs - Number of poles desired? (default: 0 --> compute all)

function [l,V] = compute_resonances(K0,K1,K2,neigs)

if nargin == 1, neigs = 0;  end
if nargin == 2, neigs = K1; end
if nargin == 3, neigs = 0;  end

if nargin == 1 | nargin == 2
  elt = K0;
  [N,nnz] = problem_size(elt);
  [K0,K1,K2] = form_operators(elt, ...
     (neigs ~= 0) & (nnz < 0.2 * N^2) & (N > 100));
end

N = length(K0);
if neigs == 0 | issparse(K0) ~= 0
  Z = zeros(N);
  I = eye(N);
else
  Z = spalloc(N,N,0);
  I = speye(N);
end
A = [K0, Z; Z, I];
B = [-K1, -K2; I, Z];

if neigs == 0
  if nargout == 1
    l = eig(A,B)*1i;
  elseif nargout == 2
    [V,D] = eig(A,B);
    l = diag(D)*1i;
  end
elseif issparse(K0)
  [L,U,P,Q] = lu(A);
  opts.isreal = 1;
  opts.disp = 0;
  if nargout == 1
    l = eigs(@cr_lufun1, 2*N, neigs, 0, opts, L, U, P, Q, B)*1i;
  elseif nargout == 2
    [V,D] = eigs(@cr_lufun1, 2*N, neigs, 0, opts, L, U, P, Q, B);
    l = diag(D)*1i;
  end
else
  [L,U] = lu(A);
  opts.isreal = 1;
  opts.disp = 0;
  if nargout == 1
    l = eigs(@cr_lufun2, 2*N, neigs, 0, opts, L, U, B)*1i;
  elseif nargout == 2
    [V,D] = eigs(@cr_lufun1, 2*N, neigs, 0, opts, L, U, P, Q, B);
    l = diag(D)*1i;
  end
end

function y = cr_lufun1(x,L,U,P,Q,B)
y = Q*(U\(L\(P*(B*x))));

function y = cr_lufun2(x,L,U,B)
y = U\(L\(B*x));
% [ll,diff] = compare_eigs(ll1, ll2, tol)
%
% Compare the results of two eigenvalue vectors (ll1 and ll2).
% Inputs:
%   ll1, ll2 - eigenvalue lists
%   tol      - matching tolerance
% Outputs:
%   ll   - all elements of ll1 within tol of some element of ll2
%   diff - diff(i) is min(abs(ll1(i) - ll2));

function [ll,diff] = compare_eigs(ll1, ll2, tol)

if size(ll1,1) == 1, ll1 = ll1.'; end
if size(ll2,1) ~= 1, ll2 = ll2.'; end
e1 = 0*ll1 + 1;
e2 = 0*ll2 + 1;
diff = min(abs(ll1*e2 - e1*ll2), [], 2);
ll   = ll1(find(diff < tol));
% plot_potential1(VV,xx)
%
% Plot the potential with continuous lines

function plot_potential1(VV,xx)

VV = [0,VV,0];
xx = [2*xx(1)-xx(2), xx, 2*xx(end)-xx(end-1)];
yy = zeros(2*length(xx),1);
WW = zeros(2*length(xx),1);
yy(1:2:end-1) = xx;
yy(2:2:end)   = xx;
WW(2:2:end-2) = VV;
WW(3:2:end-1) = VV;

plot(yy,WW,'r','LineWidth',3)
ylim([min(VV)-(max(VV)-min(VV))/3, max(VV)+(max(VV)-min(VV))/3])
% plot_potential(elt)
%
% Plot the potential

function plot_potential(elt)

nelt = length(elt);
u = zeros(problem_size(elt),1);
base = 1;
for j = 1:nelt
  order = elt(j).order;
  x = cos(pi*(order:-1:0)/order)';
  xelt = elt(j).a * (1-x)/2 + elt(j).b * (1+x)/2;
  u(base:base+order) = eval_potential(elt(j), xelt);
  base = base + order + 1;
end
plot_fields(elt,u);
% plot_fields(elt,u)
%
% Plot some function defined on the spectral element mesh

function plot_fields(elt,u)

rflag = isreal(u);
hflag = ishold;
N = problem_size(elt);

hold on;
cla;
base = 1;
for j = 1:length(elt)
  order = elt(j).order;
  x = cos(pi*(order:-1:0)/order)';

  % Plot points  
  xelt = elt(j).a * (1-x)/2 + elt(j).b * (1+x)/2;
  uelt = u(base:base+order);
  plot(xelt, real(uelt), 'r.', 'markersize', 16);
  if ~rflag, plot(xelt, imag(uelt), 'b.', 'markersize', 16); end
  
  % Plot interpolant
  xxelt = linspace(elt(j).a, elt(j).b);
  [pc,s,mu] = polyfit(xelt,uelt,order);
  uuelt = polyval(pc, (xxelt-mu(1))/mu(2));
  hr = line(xxelt, real(uuelt));
  set(hr, 'Color', 'r');
  if ~rflag
    hi = line(xxelt, imag(uuelt));
    set(hi, 'Color', 'b'); 
  end

  if length(u) > N
    base = base + order + 1;
  else
    base = base + order;
  end
end
if ~ishold, hold off; end
% [K0,K1,K2] = form_operators(elt,is_sparse)
%
% Form pseudospectral discretization of
%   ((-D^2 + V) + z^2) u = 0  on  (a,b)
%   (D-z) u = 0               at  a
%   (D+z) u = 0               at  b
% The matrices K0, K1, K2 are such that
%   K(l) = K0 + z*K1 + z^2*K2
% is the discretized operator.

function [K0,K1,K2] = form_operators(elt,is_sparse)

if nargin < 2, is_sparse = 0; end

if ~isfield(elt, 'V')
  [K0,K1,K2] = form_operators_sys(elt, is_sparse);
  return
end

% Useful constants and memory allocations
nelt  = length(elt);
if is_sparse
  [N,nz] = problem_size(elt);
  K0 = spalloc(N,N,nz);
  K1 = spalloc(N,N,2);
  K2 = spalloc(N,N,N);
else
  N  = problem_size(elt);
  K0 = zeros(N);
  K1 = zeros(N);
  K2 = zeros(N);
end

% Operator on subdomains
base = 1;
for j = 1:nelt
  
  % Element contribution (including matching conditions)
  order   = elt(j).order;
  [D,x]   = cheb(order);
  lelt(j) = elt(j).b - elt(j).a;
  xelt    = elt(j).a * (1-x)/2 + elt(j).b * (1+x)/2;
  K0elt   = -4*D^2/lelt(j)^2 + diag(eval_potential(elt(j), xelt));
  K0elt(1,:)   = -2*D(1,:)/lelt(j);
  K0elt(end,:) =  2*D(end,:)/lelt(j);
  
  % Assembly
  I             = base:base+order;
  Iint          = base+1:base+order-1;
  K0(I,I)       = K0(I,I) + K0elt;
  K2(Iint,Iint) = eye(order-1);
  
  base = base + order;
end

% Radiation boundary conditions
K1(1,1) = 1;
K1(N,N) = 1;


%
% Compute D = differentiation matrix, x = Chebyshev grid.  Adapted
% from the routine in Trefethen's "Spectral Elements in MATLAB"
%
function [D,x] = cheb(N)

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
x = x(N+1:-1:1);
D = D(N+1:-1:1,N+1:-1:1);
% N = problem_size(elt)
%
% Get total number of points in the mesh

function [N,nz] = problem_size(elt)

N = 1;
for k = 1:length(elt)
  N = N + elt(k).order;
end

if nargout == 2
  nz = 0;
  for k = 1:length(elt)
    nz = nz + (elt(k).order+1)^2;
  end
end
% [Vx] = eval_potential(eltj, x)
%
% Evaluate the potential associated with the mesh element eltj at x.

function [Vx] = eval_potential(eltj, x)

if strcmp(eltj.Vtype, 'function')
  Vx = feval(eltj.V, x, eltj.args{:});
elseif strcmp(eltj.Vtype, 'spline')
  Vx = ppval(eltj.V, x);
elseif strcmp(eltj.Vtype, 'constant')
  Vx = eltj.V + 0*x;
else
  error('Unknown potential type');
end

