% the programm is designed for figure 2.8,
% fig28.m
% **************************************************************************
% This program implements equation (2.26) using the waterfilling algorithm.
% The equation has been implemented in its determinant form as it is more convenient.
% output=capacity_water(SNR,M)
% SNR is signal-to-noise ratio; 
% M -> number of antennas for a M x M system
% **************************************************************************
% Monte-Carlo runs
       clear
       M = 4;
 testnum = 2^(15);  
     SNR = 0 : 20;
    Lsnr = length(SNR);
    for nl = 1 : Lsnr;
        snr = 10^(SNR(nl)/10);
        for K=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(K)=log2(a);   
        end
        a_var(nl) = sum(y)/testnum;
    end
%   save capacity.mat bvarible -ascii; 
    load capacity.mat -ascii,
    b_var = capacity;
    figure(1);
    plot(SNR,b_var,'k-p');
    hold on;
    plot(SNR,a_var,'k-o');
    xlabel('SNR(dB)');
    ylabel('Ergodic capacity(bit/s/Hz)');
    legend('channel unknown to Tx','channel known to Tx');
    hold off