clear
clc
M=4;
N=1;L=16;
xx=401;yy=401;
Pt=1;
snr_dB=20;
snr=10^(snr_dB/10);
I=eye(M,M);
dota=26;
sigma=8;
Pout=10/100;
for x=1:dota:801
    for y=1:dota:801
        d=sqrt((x-xx).^2+(y-yy).^2);
        for count=1:1000
            H=[];
            for i=1:N
                u(i)=10*log10(Pt/(N*L))-4*10*log10(d/1000);
                s_dB(i)=u(i)+sigma*randn(1,1);
                s(i)=10^(s_dB(i)/10);
                Hi=sqrt(s(i)/2)*randn(M,L)+sqrt(-1)*sqrt(s(i)/2)*randn(M,L);
                H=[H Hi];
            end
            C(count)=log2(det(I+(snr/(N*L))*H*H'));
        end
        C=abs(C);
        C=sort(C);
        Cout(floor(x/dota)+1,floor(y/dota)+1)=C(Pout*count+1);
    end
end
[m,n]=size(Cout);
x=1:dota:m*dota;
y=1:dota:n*dota;
mesh(x-401,y-401,Cout)
CAS_min=min(min(Cout));
CAS_max=max(max(Cout));
axis([-400 400 -400 400 0 CAS_max])
Cout_mean=sum(sum(Cout))/(m*n)

xlabel('x')
ylabel('y')
zlabel('outage capacity of CAS')