clc
clear all

B=15;
BER=10^(-4);
K=1.5/(logm(0.2/BER));
SNR=[8,10];
k=log2(1+K.*SNR);
r1=10;
r2=12;
x=0;
y=1;
z=1;
% b=[0,0];
% c=x+y*(sum(b));
% p=r.*k.*b-b.*c;

a(1)=0.001;
for i=1:300
    a(2)=0.001;
    for j=1:300
        eigenvalue(1)=1-10/3*(2*a(1)*k(1)-a(1)*k(2)+2*a(2)*k(2)-a(2)*k(1))+10/3*((a(1)*(k(2)-2*k(1))*a(2)*(k(1)-2*k(2))+(2*a(2)*k(2)-a(2)*k(1)-2*a(1)*k(1)+a(1)*k(2))^2)^(0.5));
        eigenvalue(2)=1-10/3*(2*a(1)*k(1)-a(1)*k(2)+2*a(2)*k(2)-a(2)*k(1))-10/3*((a(1)*(k(2)-2*k(1))*a(2)*(k(1)-2*k(2))+(2*a(2)*k(2)-a(2)*k(1)-2*a(1)*k(1)+a(1)*k(2))^2)^(0.5));
        if (eigenvalue(1)<1)&&(eigenvalue(1)>-1)&&(eigenvalue(2)<1)&&(eigenvalue(2)>-1)
            a(2)=a(2)+0.002;
        end
    end
    if a(2)>0
        a(2)=a(2)-0.002;
    else
        a(2)=0;
    end
    f2(i)=a(2);
    f1(i)=a(1);
    a(1)=a(1)+0.002;
end

figure(1)
plot(f1,f2,'--r')
axis([0,0.25,0,0.25])
hold on

B=15;
BER=10^(-4);
K=1.5/(logm(0.2/BER));
SNR=[10,10];
k=log2(1+K.*SNR);
r1=10;
r2=12;
x=0;
y=1;
z=1;
% b=[0,0];
% c=x+y*(sum(b));
% p=r.*k.*b-b.*c;

a(1)=0.001;
for i=1:300
    a(2)=0.001;
    for j=1:300
        eigenvalue(1)=1-10/3*(2*a(1)*k(1)-a(1)*k(2)+2*a(2)*k(2)-a(2)*k(1))+10/3*((a(1)*(k(2)-2*k(1))*a(2)*(k(1)-2*k(2))+(2*a(2)*k(2)-a(2)*k(1)-2*a(1)*k(1)+a(1)*k(2))^2)^(0.5));
        eigenvalue(2)=1-10/3*(2*a(1)*k(1)-a(1)*k(2)+2*a(2)*k(2)-a(2)*k(1))-10/3*((a(1)*(k(2)-2*k(1))*a(2)*(k(1)-2*k(2))+(2*a(2)*k(2)-a(2)*k(1)-2*a(1)*k(1)+a(1)*k(2))^2)^(0.5));
        if (eigenvalue(1)<1)&&(eigenvalue(1)>-1)&&(eigenvalue(2)<1)&&(eigenvalue(2)>-1)
            a(2)=a(2)+0.002;
        end
    end
    if a(2)>0
        a(2)=a(2)-0.002;
    else
        a(2)=0;
    end
    f2(i)=a(2);
    f1(i)=a(1);
    a(1)=a(1)+0.002;
end

figure(1)
plot(f1,f2,'-.*b')
axis([0,0.25,0,0.25])

B=15;
BER=10^(-4);
K=1.5/(logm(0.2/BER));
SNR=[10,8];
k=log2(1+K.*SNR);
r1=10;
r2=12;
x=0;
y=1;
z=1;
% b=[0,0];
% c=x+y*(sum(b));
% p=r.*k.*b-b.*c;

a(1)=0.001;
for i=1:300
    a(2)=0.001;
    for j=1:300
        eigenvalue(1)=1-10/3*(2*a(1)*k(1)-a(1)*k(2)+2*a(2)*k(2)-a(2)*k(1))+10/3*((a(1)*(k(2)-2*k(1))*a(2)*(k(1)-2*k(2))+(2*a(2)*k(2)-a(2)*k(1)-2*a(1)*k(1)+a(1)*k(2))^2)^(0.5));
        eigenvalue(2)=1-10/3*(2*a(1)*k(1)-a(1)*k(2)+2*a(2)*k(2)-a(2)*k(1))-10/3*((a(1)*(k(2)-2*k(1))*a(2)*(k(1)-2*k(2))+(2*a(2)*k(2)-a(2)*k(1)-2*a(1)*k(1)+a(1)*k(2))^2)^(0.5));
        if (eigenvalue(1)<1)&&(eigenvalue(1)>-1)&&(eigenvalue(2)<1)&&(eigenvalue(2)>-1)
            a(2)=a(2)+0.002;
        end
    end
    if a(2)>0
        a(2)=a(2)-0.002;
    else
        a(2)=0;
    end
    f2(i)=a(2);
    f1(i)=a(1);
    a(1)=a(1)+0.002;
end

figure(1)
plot(f1,f2,'-g')
axis([0,0.25,0,0.25])


B=15;
BER=10^(-4);
K=1.5/(logm(0.2/BER));
SNR=[7,7];
k=log2(1+K.*SNR);
r1=10;
r2=12;
x=0;
y=1;
z=1;
% b=[0,0];
% c=x+y*(sum(b));
% p=r.*k.*b-b.*c;

a(1)=0.001;
for i=1:300
    a(2)=0.001;
    for j=1:300
        eigenvalue(1)=1-10/3*(2*a(1)*k(1)-a(1)*k(2)+2*a(2)*k(2)-a(2)*k(1))+10/3*((a(1)*(k(2)-2*k(1))*a(2)*(k(1)-2*k(2))+(2*a(2)*k(2)-a(2)*k(1)-2*a(1)*k(1)+a(1)*k(2))^2)^(0.5));
        eigenvalue(2)=1-10/3*(2*a(1)*k(1)-a(1)*k(2)+2*a(2)*k(2)-a(2)*k(1))-10/3*((a(1)*(k(2)-2*k(1))*a(2)*(k(1)-2*k(2))+(2*a(2)*k(2)-a(2)*k(1)-2*a(1)*k(1)+a(1)*k(2))^2)^(0.5));
        if (eigenvalue(1)<1)&&(eigenvalue(1)>-1)&&(eigenvalue(2)<1)&&(eigenvalue(2)>-1)
            a(2)=a(2)+0.002;
        end
    end
    if a(2)>0
        a(2)=a(2)-0.002;
    else
        a(2)=0;
    end
    f2(i)=a(2);
    f1(i)=a(1);
    a(1)=a(1)+0.002;
end

figure(1)
plot(f1,f2,'-.y')
axis([0,0.25,0,0.25])

B=15;
BER=10^(-4);
K=1.5/(logm(0.2/BER));
SNR=[8,8];
k=log2(1+K.*SNR);
r1=10;
r2=12;
x=0;
y=1;
z=1;
% b=[0,0];
% c=x+y*(sum(b));
% p=r.*k.*b-b.*c;

a(1)=0.001;
for i=1:300
    a(2)=0.001;
    for j=1:300
        eigenvalue(1)=1-10/3*(2*a(1)*k(1)-a(1)*k(2)+2*a(2)*k(2)-a(2)*k(1))+10/3*((a(1)*(k(2)-2*k(1))*a(2)*(k(1)-2*k(2))+(2*a(2)*k(2)-a(2)*k(1)-2*a(1)*k(1)+a(1)*k(2))^2)^(0.5));
        eigenvalue(2)=1-10/3*(2*a(1)*k(1)-a(1)*k(2)+2*a(2)*k(2)-a(2)*k(1))-10/3*((a(1)*(k(2)-2*k(1))*a(2)*(k(1)-2*k(2))+(2*a(2)*k(2)-a(2)*k(1)-2*a(1)*k(1)+a(1)*k(2))^2)^(0.5));
        if (eigenvalue(1)<1)&&(eigenvalue(1)>-1)&&(eigenvalue(2)<1)&&(eigenvalue(2)>-1)
            a(2)=a(2)+0.002;
        end
    end
    if a(2)>0
        a(2)=a(2)-0.002;
    else
        a(2)=0;
    end
    f2(i)=a(2);
    f1(i)=a(1);
    a(1)=a(1)+0.002;
end

figure(1)
plot(f1,f2,'om')
axis([0,0.25,0,0.25])
legend('r1=8dB,r2=10dB','r1=10dB,r2=10dB','r1=10dB,r2=8dB','r1=7dB,r2=7dB','r1=8dB,r2=8dB')
xlabel('a1')
ylabel('a2')