function [mbest] = eigenselect(eigen);
%BREAKPOINT  
% 
%    EIGENSELECT(series) finds optimal choice of # of eigenvalues
%    to retain based on estimate of noise background trend to
%    log eigenvalue spectrum
%
%    restrict to space spanned by first 0.5*M eigenvalues (to avoid
%    non-linear regime at high eigenvalue number)
%
     N=length(eigen);
     M=round(N/2);
%    establish backround slope using first M eigenvalues
     series=log(eigen(1:M));
     N=length(series);
     Y0=series;
     xx0=(1:N)';
     wons0=ones(N,1);
     X0 = [xx0 wons0];
    [B0,BINT0,R,RINT,STATS] = regress(Y0,X0);% don't skew the fit with 1 and2
     trend=B0(1)*xx0+B0(2);
     tol=1.0*std(R);
%    Keep the set of first K eigenvalues where Kth eigenvalue
%    is the first to drop below background trend within indicated
%    tolerance
%
%    use tolerance of 1 standard deviation of trend residual
%   
     mbest=1;
     flag=0;
     for m=1:N
     if (series(m)>trend(m)+tol)
         if (flag==0) 
             mbest=m+1;
         else
         end
     else
         flag=1;
     end
     end

         
