M = input('M = ');
T = input('T = ');
B_min = input('B_min =');
B_max = input('B_max =');
Bs = B_min:B_max;
k = 0;

for index = 1:length(Bs);
    B = Bs(index);
    B_(index) = mod(B,M);
    b(index) = (B-B_(index))/M;
    
    if B_(index) <= b(index)*M/(T+b(index));
        disp('R = T/T+b')
        R(index) = T/(T+b(index));
        x(index) = T+b(index);
        u(index) = M*b(index);
        n(index) = M*x(index) - 2*u(index);
        R_Macro(index) = (u(index)+n(index))/(M*x(index));
    else
        if T ~= b(index);
            disp('R = PEC')
            R(index) = (M*(T+b(index)+1) - B)/(M*(T+b(index)+1));
            x(index) = T+b(index)+1;
            u(index) = B;
            n(index) = M*x(index) - 2*u(index);
            R_Macro(index) = (u(index)+n(index))/(M*x(index));
        else
            disp('T = Tmin = b')
            R(index) = (M-B_(index))/M;
            x(index) = 1;
            u(index) = M-B_(index);
            n(index) = 0;
            R_Macro(index) = (u(index)+n(index))/(M*x(index));
        end
    end
    
    %{
    if B_(index) <= M/2;
        disp('B_ <= M/2')
         if T <= b(index)*(M-B_(index))/B_(index);
             disp('T small, R = T/T+b')
             R(index) = T/(T+b(index));
             x(index) = T+b(index);
             u(index) = M*b(index);
             n(index) = M*x(index) - 2*u(index);
             R_Macro(index) = (u(index)+n(index))/(2*u(index)+n(index)+k*M);
         else
             disp('R = PEC')
             R(index) = (M*(T+b(index)+1) - B)/(M*(T+b(index)+1));
             x(index) = T+b(index)+1;
             u(index) = B;
             n(index) = M*x(index) - 2*u(index);
             R_Macro(index) = (u(index)+n(index))/(2*u(index)+n(index)+k*M);
         end
    else
        disp('B_ > M/2')
        if T == b(index);
            disp('T min, R = M-B_/M')
            R(index) = (M-B_(index))/M;
            x(index) = 1;
            u(index) = M-B_(index);
            n(index) = M*x(index) - 2*u(index);
            R_Macro(index) = (u(index)+n(index))/(2*u(index)+n(index)+k*M);
        else
            disp('R = PEC')
            R(index) = (M*(T+b(index)+1) - B)/(M*(T+b(index)+1));
            x(index) = T+b(index)+1;
            u(index) = B;
            n(index) = M*x(index) - 2*u(index);
            R_Macro(index) = (u(index)+n(index))/(2*u(index)+n(index)+k*M);
        end
    end
    %}
end

R_SCo = M*T./(M*T + Bs);
R_SCo(Bs > M*T) = NaN;
R_RLC = (M*(T+1)-Bs)./(M*(T+1));

%% Plotting
figure
LineWidth = 4;
FontSize = 24;
FontName = 'Times New Roman';
set(gca,'FontSize',FontSize-4);
plot(Bs,R,'s-.','LineWidth',LineWidth)
hold on
plot(Bs,R_SCo,'ro-.','LineWidth',LineWidth)
plot(Bs,R_RLC,'k+','LineWidth',LineWidth)
set(0,'DefaultLineMarkerSize',10);
grid on
title(['M = ' num2str(M) ', T = ' num2str(T)],'FontSize',FontSize,'FontName',FontName)
xlabel('Burst Length (B)','FontSize',FontSize,'FontName',FontName)
ylabel('Rate (R)','FontSize',FontSize,'FontName',FontName)
legend('Optimal Codes','SCo Codes','MDP Codes')