clear all;


N=450;% Size of Hamiltonian
samples=30;% Number of runs per iteration. "samples" number of znorm renditions will be written to output file
runs=100000;% Number of iterations. Final size of output file will contain "runs*samples" number of znorm entries. 
% Setting the value of damping 
Ktwidle=1.5;
eigen=zeros(N,samples); % Initializing the eigen matrix
    
mask1=zeros(N);% Creating hamiltonian mask for diagonal terms
mask2=zeros(N);% Creating hamiltonian mask for off-diagonal terms
            for i=1:N,
                for j=1:N,
                    if(i==j)
                        mask1(i,j)=1/2;
                    elseif(i<j)
                        mask2(i,j)=1;
                    elseif (i>j)
                        0;
                    end;
                end;
            end;


for p=1:runs,
        p
        for cnt=1:samples,
            hamilt=mask1.*normrnd(0,1,N,N)+mask2.*normrnd(0,sqrt(0.5),N,N);
            hamilt=hamilt+transpose(hamilt);% creating the resultant hamiltonian
            eigen(:,cnt)=eig(hamilt);% computing eigenvalues of hamiltonian
        end;
        
        Elevel=(N/(2*pi))*(pi+2*asin(eigen./sqrt(2*N))+2.*(eigen./sqrt(2*N)).*sqrt(2*N-eigen.^2)/sqrt(2*N))-N/2;% Mapping the eigenvalues from the semi-circle to have uniform spacing
        
        Elevel(1,:)=-9.99*ones(1,samples); % removing the outliers 
        Elevel(N,:)=9.99*ones(1,samples); % removing the outliers

        green=(1i/pi)./(Elevel+1i*Ktwidle);
        CoupleMatrix1=normrnd(0,1,N,samples);% setting up couple matrix of port 1
        CoupleMatrix2=normrnd(0,1,N,samples);% setting up couple matrix of port 2
        
        z11=sum(CoupleMatrix1.*green.*CoupleMatrix1);
        z12=sum(CoupleMatrix1.*green.*CoupleMatrix2);
        z21=sum(CoupleMatrix2.*green.*CoupleMatrix1);
        z22=sum(CoupleMatrix2.*green.*CoupleMatrix2);
        
        fid=fopen('RMTZ_1_5loss.txt','a+');
        result=[real(z11);imag(z11);real(z12);imag(z12);real(z21);imag(z21);real(z22);imag(z22)];
        fprintf(fid,'%8.5f\t %8.5f\t %8.5f\t %8.5f\t %8.5f\t %8.5f\t %8.5f\t %8.5f\t \n',result);
        fclose(fid);
end;