function cc=invnrp(A,mu);

%                                                     
% INVNRP estimates (if possible), for a given square 
%        matrix  A  and an interior point  mu  of  
%        the numerical range  F(A),  a unit vector  
%        VECTOR  such that  VECTOR* A VECTOR = mu.
%
% Authors: Ch. Chorianopoulos and P. Psarrakos,         
%          Department of Mathematics,                   
%          National Technical University of Athens.   
%
%  DATE:   December 4, 2009.                             
%                                                       

n=max(size(A)); i=sqrt(-1); B=A-mu*eye(n,n); 

tic

%%%%%%%%%%%%%%%%%%%  IMAGINARY  POINT  ib  %%%%%%%%%%%%%%

H=(B+B')/2;
[X,D]=eig(H);
[s,w] = sort(real(diag(D)));
lamin=s(1);
lamax=s(n);
xmin=X(:,w(1))/norm(X(:,w(1)),2); 
xmax=X(:,w(n))/norm(X(:,w(n)),2); 

if abs(lamin+lamax)>0.000001
   tx=(lamin+sqrt(abs(lamin*lamax)))/(lamin+lamax);
else
   tx=1/2;
end

xb=((1-tx)*xmin+tx*xmax)/norm((1-tx)*xmin+tx*xmax,2);
b=xb'*B*xb;

if abs(b)<0.0000001
    b=imag(b)
    VECTOR=xb
    ESTIMATION=VECTOR'*A*VECTOR
    absolute_error=abs(ESTIMATION-mu)
end

%%%%%%%%%%%%%%%%  REAL  POINT  a  %%%%%%%%%%%%%%%%%%%%%%%

K=(B-B')/(2*i);
[Y,D]=eig(K);
[s,w]=sort(real(diag(D)));
lamin=s(1);
lamax=s(n);
ymin=Y(:,w(1))/norm(Y(:,w(1)),2);
ymax=Y(:,w(n))/norm(Y(:,w(n)),2); 

if abs(lamin+lamax)>0.000001
   ty=(lamin+sqrt(abs(lamin*lamax)))/(lamin+lamax);
else
   ty=1/2;
end

xa=((1-ty)*ymin+ty*ymax)/norm((1-ty)*ymin+ty*ymax,2);
a=xa'*B*xa;

if abs(b)>0.0000001 & abs(a)<0.0000001
    a=real(a)
    VECTOR=xa
    ESTIMATION=VECTOR'*A*VECTOR
    absolute_error=abs(ESTIMATION-mu)
end

%%%%%%%%%%%%%%%%  BOUNDARY  POINT  q  %%%%%%%%%%%%%%%%%%%%%%%

if abs(a)>0.0000001 & abs(b)>0.0000001   

   fi=pi+angle(a+b*abs(a/b));
   z=exp(i*(-fi));                 % the angular of the rotation
   B1=z*B;                         % the rotation of matrix 
   B2=(B1+B1')/2;                  % the hermitian part of the rotation 
   [X,D]=eig(B2);                  % eigenvectors and eigenvalues
   [s,w] = sort(real(diag(D)));
   xq=X(:,w(n))/norm(X(:,w(n)),2); % unitary eigenvector of the maximum eigenvalue
   q=xq'*B*xq;                     % boundary point
   
k=1;
fii=fi;

while abs(exp(i*fi)-exp(i*angle(q)))>2*sin(pi/8) & k<20
    k=k+1;
    fii=fii+(sign(angle(exp(i*fi)/exp(i*angle(q))))*pi)/(2^(k+1));
    z=exp(i*(-fii));                % the angular of the rotation
    B1=z*B;                         % the rotation of the matrix 
    B2=(B1+B1')/2;                  % the hermitian part of the rotation 
    [X,D]=eig(B2);                  % eigenvectors and eigenvalues
    [s,w] = sort(real(diag(D)));
    xq=X(:,w(n))/norm(X(:,w(n)),2); % unitary eigenvector of the maximum eigenvalue
    q=xq'*B*xq;                     % boundary point   
end

c=real(imag(b)*real(q))/(imag(b)-imag(q));
gwnia=angle(q-b);
G=(B-c*eye(n,n))*exp(-i*gwnia); 
pphi0=angle(xq'*G*xb-xb.'*G.''*xq.'');
t0=real( ( sqrt(((exp(i*pphi0)*xb'*G*xq+exp(-i*pphi0)*xq'*G*xb))^2-4*(-abs(b-c))*(abs(q-c))))-(exp(i*pphi0)*xb'*G*xq+exp(-i*pphi0)*xq'*G*xb) ) / (2*abs(q-c)) ;
xc=exp(-i*pphi0)*xb+t0*xq;
xc=xc/norm(xc,2);

phi0=angle(xc'*B*xa-xa.'*B.''*xc.'');
t0=real( (  sqrt(((exp(i*phi0)*xa'*B*xc+exp(-i*phi0)*xc'*B*xa))^2-4*a*c)-(exp(i*phi0)*xa'*B*xc+exp(-i*phi0)*xc'*B*xa) ) / (2*c) );

VECTOR=exp(-i*phi0)*xa+t0*xc;
VECTOR=VECTOR/norm(VECTOR,2)
ESTIMATION=VECTOR'*A*VECTOR
absolute_error=abs(ESTIMATION-mu)

toc

number_of_iterations_for_point_q = k

if k==20
    disp('The point is near the boundary OR in the exterior of numerical range, OR the boundary has flat portions!!')
end

end
