% MATHLAB/OCTAVE HINTS
% PHYSICS 6562
% STATISTICAL MECHANICS

% Matlab example for random matrix exercise
% This hints file covers only the computational parts (a & f) of the exercise:
% Don't neglect the analytical portions!

% GENERATING RANDOM GOE MATRICES

% Generating M members of the GOE ensemble of NxN matrices, and 
% calculating the difference between the eigenvalues 
% in the middle of the range

% Increase M, N, and Nbins as appropriate

M = 20;
N = 4;
Nbins = 8;

% randn gives Gaussian random numbers; randn(N) generates a NxN matrix
% Mat' is the transpose

differences = zeros(M,1);	% Initialize array of differences
M11 = zeros(M,1);	% Initialize array of 11 entries

for m = 1:M
    Mat = randn(N);
    Ms = Mat + Mat';
    lambda = sort(eig(Ms));
    differences(m) = lambda(N/2+1)-lambda(N/2);
    M11(m) = Ms(1,1);
end

% Divide out by mean value of the splittings
diffAve = mean(differences)

% hist generates a histogram of differences with Nbins bins. 
% Change Nbins to 50 or so when you do production runs. 

% Octave allows one to use the third argument to normalize the area of
% the histogram (say, to one).
% Matlab doesn't normalizes the histograms 

hist(differences/diffAve,Nbins)

% PLOTTING THEORY OVER HISTOGRAM

% Demonstrated with histogram for diagonal element
% First plot histogram for M11
hist(M11,Nbins)

% "hold on" and "hold off" allow one to combine curves
hold on

% Now, plot expected Gaussian fit

% define x curve
x = min(M11)-1:0.01:max(M11)+1;
M11Range = max(M11)-min(M11)

% y(x) is Gaussian of RMS width sigma = 2 (diagonal element doubled in size)
% expected width of diagonal element
sigma11 = 2;
y = (1/sqrt(2*pi*sigma11^2))*exp(-x.*x/(2*sigma11^2));

% Histogram multiplies height by number of entries, divides it by Nbins
normalization = M*M11Range/Nbins
plot(x,normalization*y)

% Reset graphics so next curve is not overlayed on these two
hold off

% GENERATING RANDOM +-1 MATRICES

% Symmetric matrix with integer values +-1 with 50/50 probability
% rand gives flat distributions of random numbers between 0 and 1
% generate random matrix of +- 1
MatPM = sign(2*rand(N) -1)

% Allocate space for NxN matrix Ms
MsPM = zeros(N);

% Symmetrize: copy top half of the matrix (diagonal copied twice)
for i=1:N
  for j=i:N
    MsPM(i,j) = MatPM(i,j);
    MsPM(j,i) = MatPM(i,j);
  end
end  

% Show the matrix: check that it worked
MsPM