% the programm is designed for figure 2.10,
% fig210.m
% Monte-Carlo runs
       clear
       M = 4;
 testnum = 2^(5);
 probnum = 10^3;
  outnum = 10^2;
     SNR = 0 : 20;
    Lsnr = length(SNR);
    for nl = 1 : Lsnr;
        snr = 10^(SNR(nl)/10);
        for lp = 1 : probnum;
            for kk = 1 : testnum,
                T = randn(M,M)+j*randn(M,M);
                T = T/sqrt(2);
                I = eye(M);
                eigen = eig(T*T'); % extract eigenvalues
                gamma = zeros(M,1);
                gamma = watfill(snr,T,M);
                a=I+(snr/M)*diag(gamma).*diag(eigen);% solve eqn(2.26)using the determinant form
                a = det(a);
	            y(kk) = log2(a);
            end
            a_var(lp) = sum(y)/testnum;
        end
        d_var = sort(a_var);
        D_var(nl) = d_var(outnum); 
    end
%   save capacity29.mat bvar -ascii; 
    load capacity29.mat -ascii,
    b_var = capacity29;
    figure(1);
    plot(SNR,b_var,'k-p');
    hold on;
    plot(SNR,D_var,'k-o');
    xlabel('SNR(dB)');
    ylabel('10% outage capacity(bit/s/Hz)');
    legend('channel unknown to Tx','channel known to Tx');
    hold off