%MATLAB code to create Figure 7(b) in our IEEE Globecom 2008 paper:
%   J.Katto, K.Ogura, Y.Akae, T.Fujikawa, K.Kaneko, S.Zhou: 
%   "Simple Model Analysis and Performance Tuning of Hybrid TCP Congestion 
%   Control", IEEE Globecom 2008, Dec.2008.
%
%   Just type 'GCOM002'. 
%   Default parameter settings are
%      Bandwidth: 100Mbps
%      RTT: 40msec
%      Packet size: 1500bytes
%   You can modify these parameters and can evaluate several TCP performances.
%
%   Copyright (c) 2008, Katto Lab., Waseda University.
%   All rights reserved.
%
clear all;

% Two Flow Model

%% Parameter definition
%%% X-axis's Scale
MAX = 6;
%%% Packet Size
Ps = 1500 * 8;
%%% Bandwidth
BE = 100 * 10^6;
%%% RTT
RTT = 40 * 10^(-3);
%%% BDP
W = BE * RTT / Ps;
%%% Border Probability 
P1 = 1 / (3/8 * W^2);   %% w=W
P2 = 1 / (3/32 * W^2);  %% w=W/2

%%% Vector Initialization
scale = 16;
loss = zeros(1,MAX*scale);
hybrid = zeros(1,MAX*scale);
ctcp = zeros(1,MAX*scale);
areno = zeros(1,MAX*scale);
t1_ctcp = zeros(1,MAX*scale);
t2_ctcp = zeros(1,MAX*scale);
t3_ctcp = zeros(1,MAX*scale);
t1_yeah = zeros(1,MAX*scale);
t2_yeah = zeros(1,MAX*scale);
t3_yeah = zeros(1,MAX*scale);
t4_yeah = zeros(1,MAX*scale);
X = zeros(1,MAX*scale);

%% Calculation
for i=1:MAX*scale
    %%% Parameters
    %%%% PLR (packet loss rate)
    p = 10^(-(MAX-i/scale+1/scale));
    %%%% cwnd
    w = sqrt(8.0/(3.0*p));
    %%%% X-axis
    X(i) = -(MAX-i/scale+1/scale);
    %%%% Timeout penalty
    to = 2 * RTT * (1+p+2*p^2+4*p^3+8*p^4+16*p^5+32*p^6) / (1-p);
    
    %%% CTCP
    MAXTIME = 10^6;
    if( W>w ) 
        %%%% t1 for CTCP (cwnd(0)=w/2)
        w_ctcp = zeros(1,MAXTIME);
        w_ctcp_cwnd = zeros(1,MAXTIME);
        w_ctcp_dwnd = zeros(1,MAXTIME);
        w_ctcp(1) = w*0.5;
        w_ctcp_cwnd(1) = w_ctcp(1);
        w_ctcp_dwnd(1) = 0;
        for j=2:MAXTIME
            if ( w_ctcp(j-1) > (-(j-2)+(W-w/2)) ) 
                break;
            else
                w_ctcp_cwnd(j) = w_ctcp_cwnd(j-1) + 1;
                w_ctcp_dwnd(j) = w_ctcp_dwnd(j-1) + 0.125 * w_ctcp(j-1)^0.75 -1;
                if( w_ctcp_dwnd(j) < 0) 
                    w_ctcp_dwnd(j) = 0;
                end
                w_ctcp(j) = w_ctcp_dwnd(j) + w_ctcp_cwnd(j);
            end
        end
        index = j;
        t1_ctcp(i) = MAXTIME;
        for j=2:index
            if( w_ctcp(j-1) < (-(j-2)+(W-w/2)) && w_ctcp(j) >= (-(j-1)+(W-w/2)) )
                d1 = (-(j-2)+(W-w/2)) - w_ctcp(j-1);
                d2 = w_ctcp(j) - (-(j-1)+(W-w/2));
                t1_ctcp(i) = (j-1) + d1/(d1+d2) - 1;
                break;
            end
        end
        %%%% t2 for CTCP (cwnd(0)=(W-w)/2)
        w_ctcp2 = zeros(1,MAXTIME);
        w_ctcp2_cwnd = zeros(1,MAXTIME);
        w_ctcp2_dwnd = zeros(1,MAXTIME);
        w_ctcp2(1) = (W-w)/2;
        w_ctcp2_cwnd(1) = w_ctcp2(1);
        w_ctcp2_dwnd(1) = 0;
        for j=2:MAXTIME
            if ( w_ctcp2(j-1) > (-(j-2)+(W-w/2))) 
                break;
            else
                w_ctcp2_cwnd(j) = w_ctcp2_cwnd(j-1) + 1;
                w_ctcp2_dwnd(j) = w_ctcp2_dwnd(j-1) + 0.125 * w_ctcp2(j-1)^0.75 -1;
                if( w_ctcp2_dwnd(j) < 0) 
                    w_ctcp2_dwnd(j) = 0;
                end
                w_ctcp2(j) = w_ctcp2_dwnd(j) + w_ctcp2_cwnd(j);
            end
        end
        index = j;
        t2_ctcp(i) = MAXTIME;
        for j=2:index
            if( w_ctcp2(j-1) < (-(j-2)+(W-w/2)) && w_ctcp2(j) >= (-(j-1)+(W-w/2)) )
                d1 = (-(j-2)+(W-w/2)) - w_ctcp2(j-1);
                d2 = w_ctcp2(j) - (-(j-1)+(W-w/2));
                t2_ctcp(i) = (j-1) + d1/(d1+d2) - 1;
                break;
            end
        end
        %%%% t3 for CTCP (cwnd(0)=w/2)
        w_ctcp3 = zeros(1,MAXTIME);
        w_ctcp3_cwnd = zeros(1,MAXTIME);
        w_ctcp3_dwnd = zeros(1,MAXTIME);
        w_ctcp3(1) = w/2;
        w_ctcp3_cwnd(1) = w_ctcp3(1);
        w_ctcp3_dwnd(1) = 0;
        for j=2:MAXTIME
            if ( j > w/2) 
                break;
            else
                w_ctcp3_cwnd(j) = w_ctcp3_cwnd(j-1) + 1;
                w_ctcp3_dwnd(j) = w_ctcp3_dwnd(j-1) + 0.125 * w_ctcp3(j-1)^0.75 -1;
                if( w_ctcp3_dwnd(j) < 0) 
                    w_ctcp3_dwnd(j) = 0;
                end
                w_ctcp3(j) = w_ctcp3_dwnd(j) + w_ctcp3_cwnd(j);
            end
        end
        index = j;
        t3_ctcp(i) = j-1;
        for j=2:index
%            if( w_ctcp3(j-1) < w && w_ctcp3(j) >= w )
%                d1 = w - w_ctcp3(j-1);
%                d2 = w_ctcp3(j) - w;
%                t3_ctcp(i) = (j-1) + d1/(d1+d2) - 1;
%                break;
%            end
        end
    end

    %%% ARENO
    Winc_areno = BE / (10*10^6);
    t1_areno(i) = (W-w) / (1+Winc_areno);
    t2_areno(i) = W/2 / (1+Winc_areno);
    
    %%% YeAH-TCP
    MAXTIME = 10^6;
    stcp_coef = 1.01;
    if( W>w ) 
        %%%% t1 for YeAH-TCP
        w_yeah = zeros(1,MAXTIME);
        w_yeah(1) = W/2;
        for j=2:MAXTIME
            if ( w_yeah(j-1) > (-(j-2)+(W-w/2)) ) 
                break;
            else
                w_yeah(j) = stcp_coef * w_yeah(j-1);
            end
        end
        index = j;
        t1_yeah(i) = MAXTIME;
        for j=2:index
            if( w_yeah(j-1) < (-(j-2)+(W-w/2)) && w_yeah(j) >= (-(j-1)+(W-w/2)) )
                d1 = (-(j-2)+(W-w/2)) - w_yeah(j-1);
                d2 = w_yeah(j) - (-(j-1)+(W-w/2));
                t1_yeah(i) = (j-1) + d1/(d1+d2) - 1;
                break;
            end
        end
        %%%% t2 for YeAH-TCP
        w_yeah2 = zeros(1,MAXTIME);
        w_yeah2(1) = (W-w);
        for j=2:MAXTIME
            if ( w_yeah2(j-1) > (-(j-2)+(W-w/2)) ) 
                break;
            else
                w_yeah2(j) = stcp_coef * w_yeah2(j-1);
            end
        end
        index = j;
        t2_yeah(i) = MAXTIME;
        for j=2:index
            if( w_yeah2(j-1) < (-(j-2)+(W-w/2)) && w_yeah2(j) >= (-(j-1)+(W-w/2)) )
                d1 = (-(j-2)+(W-w/2)) - w_yeah2(j-1);
                d2 = w_yeah2(j) - (-(j-1)+(W-w/2));
                t2_yeah(i) = (j-1) + d1/(d1+d2) - 1;
                break;
            end
        end
    end
    if (W >= 7*w/4)        
        %%%% t3 for YeahTCP
        w_yeah3 = zeros(1,MAXTIME);
        w_yeah3(1) = w*7/8;
        for j=2:MAXTIME
            if ( w_yeah3(j-1) > (-(j-2)+(W-w/2)) ) 
                break;
            else
                w_yeah3(j) = stcp_coef * w_yeah3(j-1);
            end
        end
        index = j;
        t3_yeah(i) = MAXTIME;
        for j=2:index
            if( w_yeah3(j-1) < (-(j-2)+(W-w/2)) && w_yeah3(j) >= (-(j-1)+(W-w/2)) )
                d1 = (-(j-2)+(W-w/2)) - w_yeah3(j-1);
                d2 = w_yeah3(j) - (-(j-1)+(W-w/2));
                t3_yeah(i) = (j-1) + d1/(d1+d2) - 1;
                break;
            end
        end
        %%%% t4 for YeahTCP
        w_yeah4 = zeros(1,MAXTIME);
        w_yeah4(1) = 7*w/8;
        for j=2:MAXTIME
            if ( w_yeah4(j-1) >= w ) 
                break;
            else
                w_yeah4(j) = stcp_coef * w_yeah4(j-1);
            end
        end
        index = j;
        t4_yeah(i) = MAXTIME;
        for j=2:index
            if( w_yeah4(j-1) < w && w_yeah4(j) >= w )
                d1 = w - w_yeah4(j-1);
                d2 = w_yeah4(j) - w;
                t4_yeah(i) = (j-1) + d1/(d1+d2) - 1;
                break;
            end
        end
    end
    
    %%% TCP-Fusion
    DMIN = 1 * 10^(-3);
    Winc_fusion = BE * DMIN / Ps;
    if(W>w) 
        t1_fusion(i) = (W-w)/2 / (1+Winc_fusion);
        t2_fusion(i) = w/2 / (1+Winc_fusion);
    end
    
    %%% Parameter Initialization
    loss_num = 0;
    hybrid_num = 0;
    ctcp_num = 0;
    areno_num = 0;
    yeah_num = 0;
    fusion_num = 0;
    
    %%% Loss(Reno) & Hybrid
    if (p <= P1) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/4 * w * (3*w-2*W) * Ps/BE;
        hybrid_num = 3/8 * w^2 * Ps;
        hybrid_rtt = loss_rtt;
        hratio = hybrid_num / loss_num;
    elseif (p > P1 && p <= P2 )
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/4 * (2*w-W)^2 * Ps/BE;
        hybrid_num = 3/8 * w^2 * Ps + 1/4 * (W-w)^2 * Ps;
        hybrid_rtt = loss_rtt;
        hratio = hybrid_num / loss_num;
    else
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;
        hybrid_num = (W * w/2 - 3/8 * w^2) * Ps;
        hybrid_rtt = loss_rtt;
        hratio = hybrid_num / loss_num;
    end        
    loss(i) = loss_num / (loss_rtt + to) / 10^6;
    hybrid(i) = hybrid_num / (hybrid_rtt + hratio * to) / 10^6;

    %%% YeAH-TCP
    if (W<w) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/4 * w * (3*w-2*W) * Ps/BE;
        yeah_num = 3/8 * w^2 * Ps;
        yeah_rtt = loss_rtt;
        yratio = yeah_num / loss_num;
    elseif (w<W && W<=2*(7*w/8))
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/4 * (2*w-W)^2 * Ps/BE;

        index = fix(t1_yeah(i));
        S1 = 0;
        for j=1:index
            S1 = S1 + w_yeah(j);
        end
        S1 = S1 + w_yeah(index+1) * (t1_yeah(i)-index);
        w1 = -t1_yeah(i)+W-w/2;
        S1 = S1 + 1/2*(w1+W/2)*((W-w)/2-t1_yeah(i)) - 1/8*(W+w)*(W-w);
        
        yeah_num = (3/8*w^2 + S1) * Ps;
        yeah_rtt = loss_rtt;
        yratio = yeah_num / loss_num;
    elseif (2*(7*w/8)<W && W<=2*w)
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/4 * (2*w-W)^2 * Ps/BE;

        index = fix(t3_yeah(i));
        S1 = 0;
        for j=1:index
            S1 = S1 + w_yeah3(j);
        end
        S1 = S1 + w_yeah3(index+1) * (t3_yeah(i)-index);
        w3 = -t3_yeah(i)+W-w/2;
        S1 = S1 + 1/2*(w1+W/2)*((W-w)/2-t3_yeah(i)) - 1/8*(W+w)*(W-w);
        yeah_num = (3/8*w^2 + S1) * Ps;
        yeah_rtt = loss_rtt;
        yratio = yeah_num / loss_num;
    elseif (2*w<W && t3_yeah(i)<w/2)
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;

        index = fix(t3_yeah(i));
        S1 = 0;
        for j=1:index
            S1 = S1 + w_yeah3(j);
        end
        S1 = S1 + w_yeah3(index+1) * (t3_yeah(i)-index);
%        if (S1*Ps < loss_num) 
%            yeah_num = loss_num; 
%            yeah_rtt = RTT * (t3_yeah(i) + (loss_num/Ps - S1) / W);
%            yratio = 1;
%        else          
%            yeah_num = loss_num;
%            yeah_rtt = RTT * t3_yeah(i);
%            yratio = 1;
%        end        
        w1 = -t3_yeah(i)+W-w/2;
        S1 = S1 + 1/2*(w1+(W-w))*(w/2-t3_yeah(i));
        yeah_num = S1 * Ps;
        yeah_rtt = loss_rtt;
        yratio = yeah_num / loss_num;
    else
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt  = RTT * w/2;
        
        index = fix(t4_yeah(i));
        S2 = 0;
        for j=1:index
            S2 = S2 + w_yeah4(j);
        end
        S2 = S2 + w_yeah4(index+1) * (t4_yeah(i)-index);
        yeah_num = S2 * Ps;
        yeah_rtt = RTT * t4_yeah(i);
        yratio = yeah_num / loss_num;
    end        
    yeah(i) = yeah_num / (yeah_rtt + yratio * to) / 10^6;

    %%% CTCP
    if (W<w) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/4 * w * (3*w-2*W) * Ps/BE;
        ctcp_num = 3/8 * w^2 * Ps;
        ctcp_rtt = loss_rtt;
        cratio = ctcp_num / loss_num;
    elseif (w<W && 2*w>=W) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/4 * (2*w-W)^2 * Ps/BE;

        index = fix(t1_ctcp(i));
        S1 = 0;
        for j=1:index
            S1 = S1 + w_ctcp(j);
        end
        S1 = S1 + w_ctcp(index+1) * (t1_ctcp(i)-index);
        w1 = -t1_ctcp(i)+W-w/2;
        S1 = S1 + 1/2*(w1+W/2)*((W-w)/2-t1_ctcp(i)) - 1/8*(W+w)*(W-w);
        ctcp_num = (3/8*w^2 + S1) * Ps;
        ctcp_rtt = loss_rtt;
        cratio = ctcp_num / loss_num;
%    elseif (w<W && t2_ctcp(i)<w/2) 
    elseif (w<W && t1_ctcp(i)<w/2) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;
%        index = fix(t2_ctcp(i));
%        S1 = 0;
%        for j=1:index
%            S1 = S1 + w_ctcp2(j);
%        end
%        S1 = S1 + w_ctcp2(index+1) * (t2_ctcp(i)-index);
%        w1 = -t2_ctcp(i)+W-w/2;
%        S1 = S1 + 1/2*(w1+(W-w))*(w/2-t2_ctcp(i));
        index = fix(t1_ctcp(i));
        S1 = 0;
        for j=1:index
            S1 = S1 + w_ctcp(j);
        end
        S1 = S1 + w_ctcp(index+1) * (t1_ctcp(i)-index);
        w1 = -t1_ctcp(i)+W-w/2;
        S1 = S1 + 1/2*(w1+(W-w))*(w/2-t1_ctcp(i));
        ctcp_num = S1 * Ps;
        ctcp_rtt = loss_rtt;
        cratio = ctcp_num / loss_num;
    else
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;
%        index = fix(t1_ctcp(i));
%        S2 = 0;
%        for j=1:index
%            S2 = S2 + w_ctcp(j);
%        end
%        S2 = S2 + w_ctcp(index+1) * (t1_ctcp(i)-index);
%        ctcp_num = S2 * Ps;
%        ctcp_rtt = RTT * t1_ctcp(i);
        index = fix(t3_ctcp(i));
        S2 = 0;
        for j=1:index
            S2 = S2 + w_ctcp3(j);
        end
        S2 = S2 + w_ctcp3(index+1) * (t3_ctcp(i)-index);
        ctcp_num = S2 * Ps;
%        ctcp_rtt = RTT * t3_ctcp(i);
        ctcp_rtt = RTT * index;

        cratio = ctcp_num / loss_num;
    end;    
    ctcp(i) = ctcp_num / (ctcp_rtt + cratio * to) / 10^6;

    %%% ARENO
    if (W<w) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/4 * w * (3*w-2*W) * Ps/BE;
        areno_num = 3/8 * w^2 * Ps;
        areno_rtt = loss_rtt;
        aratio = areno_num / loss_num;
    elseif (w<W && 2*w>=W) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/4 * (2*w-W)^2 * Ps/BE;

        w1 = -t1_areno(i)+W-w/2;
        S1 = 1/2*(w/2+w1)*t1_areno(i);
        S1 = S1 + 1/2*(w1+W/2)*((W-w)/2-t1_areno(i));
        S1 = S1 - 1/8*(W+w)*(W-w);
        areno_num = (3/8*w^2 + S1) * Ps;
        areno_rtt = loss_rtt;
        aratio = areno_num / loss_num;
    elseif (w<W && t2_areno(i)<w/2) 
%    elseif (w<W && t1_areno(i)<w/2) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;
%        w1 = -t1_areno(i)+W-w/2;
%        S1 = 1/2*(w/2+w1)*t1_areno(i);
%        S1 = S1 + 1/2*(w1+(W-w))*(w/2-t1_areno(i));
        w1 = -t2_areno(i)+W-w/2;
        S1 = 1/2*((W-w)/2+w1)*t2_areno(i);
        S1 = S1 + 1/2*(w1+(W-w))*(w/2-t2_areno(i));
        areno_num = S1 * Ps;
        areno_rtt = loss_rtt;
        aratio = areno_num / loss_num;
    else
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;
%        areno_num = (3/8*(W-3/4*w)^2) * Ps;
%        areno_num = (W * w/2 - 3/8 * w^2) * Ps;
%        areno_rtt = loss_rtt;
%        aratio = areno_num / loss_num;
        areno_num = 3/8 * W^2/Winc_areno * Ps;
        areno_rtt = RTT * W/2/Winc_areno;
        aratio = areno_num / loss_num;
    end;    
    areno(i) = areno_num / (areno_rtt + aratio * to) / 10^6;

    %%% Fusion
    t2 = w / (2*(1+Winc_fusion));
    if (W < w) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt  = RTT * w/2 + 1/4 * w * (3*w-2*W) * Ps/BE;
        fusion_num = 3/8 * w^2 * Ps;
        fusion_rtt = loss_rtt;
        fratio = fusion_num / loss_num;
    elseif (w<W && 2*w>=W) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2 + 1/4 * (2*w-W)^2 * Ps/BE;
        t1 = (W-w) / (2*(1+Winc_fusion));
        y1 = Winc_fusion*t1 + W/2;
        S1 = 1/2*(W/2+y1)*t1 + 1/2*(y1+W/2)*((W-w)/2-t1) - 1/8*(W+w)*(W-w);
        fusion_num = (3/8*w^2 + S1) * Ps;
        fusion_rtt = loss_rtt;
        fratio = fusion_num / loss_num;
    elseif (2*w<W && t2<w/2) 
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt = RTT * w/2;
        y2 = Winc_fusion*t2 + (W-w);
        S2 = 1/2*((W-w)+y2)*t2 + 1/2*(y2+(W-w))*(w/2-t2);
        fusion_num = S2 * Ps;
        fusion_rtt = loss_rtt;
        fratio = fusion_num / loss_num;
    else
        loss_num = 3/8 * w^2 * Ps;
        loss_rtt  = RTT * w/2;
        fusion_num = (W * w/2 - 3/8 * w^2) * Ps;
        fusion_rtt = loss_rtt;
        fratio = fusion_num / loss_num;
    end
    fusion(i) = fusion_num / (fusion_rtt + fratio * to) / 10^6;
end;

%% Graphs
plot(X, loss, X, hybrid, X, ctcp, X, areno, X, yeah, X, fusion);
line([log10(P1),log10(P1)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
line([log10(P2),log10(P2)],[0,(BE*1.2)/10^6], 'LineStyle', ':'); 
set(gca,'XLim',[-MAX -1]);
set(gca,'YLim',[0 (BE)/10^6]);
legend('Loss(Reno)','Hybrid','CTCP','ARENO','YeAH','Fusion');
xlabel('Packet Loss Rate [10\^n]');
ylabel('Throughput [Mbps]');

%% CSV file output
z = zeros(7,96);
z(1,:)=X;
z(2,:)=loss;
z(3,:)=hybrid;
z(4,:)=ctcp;
z(5,:)=areno;
z(6,:)=yeah;
z(7,:)=fusion;
csvwrite('two_flow.txt',z');
