% Yi-Kai Liu
% 5/10/01

% N * N real symmetric matrices, BAND-DIAGONAL: 
% Elements within distance r of the main diagonal 
% are i.i.d. Uniform, between -1 and 1, with mean 0; 
% all other elements are 0.

% Implemented using regular matrix representation. 

% OK: Mehta normalizes by sqrt(2*N), but I use 2*sqrt(N) ??? 
% Apparently this is due to some notation trouble. In Porter's 
% book, p.197, the semicircle is normalized by 2*sqrt(N). 

% DONE: Exporting data
% FIXME: Examine debugging output

distribution_name = 'Uniform band-diagonal';

N = input('N * N matrices, choose N: ');
num_matrices = input('Number of matrices: ');

r = input('r = distance from main diag to edge of band, \nchoose r: ');
distribution_name = [distribution_name, ' (r=', num2str(r), ')'];

save_to_disk = input('Save output to disk? (1=yes, 0=no) ');
if save_to_disk == 1
   base_filename = input('Filename (without suffix): ', 's');
end;

A = zeros(N);                   % Random matrix
B = zeros(N, r+1);              % Used to generate A
d = zeros(N, num_matrices);     % Eigenvalues
sp = zeros(N-1, num_matrices);  % Spacings


for m = 1 : num_matrices
  
  % Generate random matrix--vectorized
  
  B = 2.0*rand(N, r+1) - 1.0;
  A = diag(B(:,1));   % Main diagonal
  for k = 2 : r+1
    A = A + diag( B(1:N-(k-1), k), k-1 );
      % Insert diagonals above main
  end;
  A = A + (triu(A,1))';   % Symmetrize
  
  % Normalized eigenvalues in interval [-1,1]
  d(:,m) = eig(A)./(0.7*sqrt(N));
    % CHANGED normalization factor
  
  sd = sort(d(:,m));   % Sorted eigenvalues

  % Compute spacings
  % for k = 1 : N-1
  %   sp(k,m) = sd(k+1) - sd(k);
  % end
  
  % Compute spacings--vectorized
  sp(:,m) = sd(2:N) - sd(1:N-1);
  
  sp(:,m) = sp(:,m) ./ mean(sp(:,m));   % Normalized spacings
                                        % with mean=1
  
end


% Save variables to disk
if save_to_disk == 1
   save([base_filename, '-data'],... 
        'distribution_name', 'save_to_disk', 'base_filename',... 
        'N', 'num_matrices', 'd', 'sp');
end;

% Plot distributions
plotdistrib;