function [asymptotic_factor] = FoV_asymptotic_factor(A,param)
% [asymptotic_factor] = FoV_asymptotic_factor(A,param)
%
% computes the asymptotic factor in a FoV (=field-of-value) analysis of GMRES
% applied to a linear system of the form Ax=b
%
% A : the matrix
% param : a vector of parameters of length 1 or 2 (optional)
% param(1): the maximum number of edges of the polygon approaching the FoV of A (10 by defaults)
%           the computation of this polygon uses the matlab function eigs for computing largest 
%           eigenvalues of the real part of e^{-it} A
% param(2) \in \{ 0,1,... \} (default 0): 
%           indicates in which figure the FoV should be plotted 
%           (=0 default figure, <0: no plot)
%
% This program only works if one has installed before the sc toolbox of Toby Driscoll, see his Webpage 
% and if one has added via the matlab command "addpath" the directory of this toolbox to the path
%
% Run FoV_demo to see some examples
%
% Bernd Beckermann, bbecker@math.univ-lille1.fr, April 13, 2005


if nargin ==1 
   param=[10,0];   
elseif nargin ==2
   if length(param)==1 
      param=[param,0];
   end;   
else
   error('FoV_asymptotic_factor takes only two arguments')
end;

% part 1: find the vertices of a polygon stored in vertices containing FoV(A)
slopes=[]; % the edges are first defined by straight lines z(\alpha)=a+b*\alpha, a,b complex
h=waitbar(0,'Computation of eigenvalues of real parts of e^{-it} A, please wait...'); 
for j=1:param(1)
   waitbar(j/param(1))
   rotation=exp(2*pi*i*j/param(1));
   lambda_max=max(real(eigs(((A/rotation)+(A/rotation)')/2,1,'LR'))); % computes the eigenvalue with largest real part
   a=rotation*lambda_max;
   b=i*rotation;
   slopes=[slopes;[a,b,lambda_max]];
end;
close(h);
slopes=[slopes;slopes(1,:)]; % repeat first edge
vertices=[];
for j=1:param(1)
   alpha=imag( (slopes(j+1,1)-slopes(j,1)) * conj(slopes(j+1,2)) );
   alpha=alpha/imag( slopes(j,2) * conj(slopes(j+1,2)) );
   vertices=[vertices,slopes(j,1)+slopes(j,2)*alpha];
end;

% part 2: check whether origin lies in the polygon, in this case error message
p=polygon(vertices);
inside=1;
for j=1:param(1)
   if imag(slopes(j,1)*conj(slopes(j,2)))>0 
      inside=0;
   end;
end;
if inside==1 % origin lies in the polygon
   if param(2)>=0
      if param(2)>0, figure(param(2));end;
      plot(p);
      hold on;
      plot(0,0,'ro');
      hold off;
      title(['Asymptotic convergence factor = 1, origin in polygon']);
   end;   
   disp('FoV WARNING: origin lies in polygon approaching FoV, try with higher order discretization');
   asymptotic_factor=1;
else % origin lies outside polygon
   % part 3: draw polygon and level lines of the exterior Riemann map
   f=extermap(p);
   asymptotic_factor=abs(evalinv(f,0));
   if param(2)>=0
      if param(2)>0, figure(param(2));end;
      plot(f);
      hold on;
      plot(0,0,'ro');
      title(['Asymptotic convergence factor ',num2str(asymptotic_factor)]);
      hold off;
   end;   
end;
