% tcpsim.m - congestion control simulation
% RMM, 9 Sep 06

aminit;					% initialize macros

congctrl_params;			% load congestion control parameters
N = 60;					% start with 60 identical sources
[weq, beq] = congctrl_equil(N, rho, c)

% Compute a perturbation for the initial condition
randn('state', 0); rand('seed', 2);
w0 = weq * ones(1, N/10) + 2 * (rand(1, N/10) - 0.5) * weq
b0 = beq/2;

% Compute the equlibrium values
[t1, y1] = ode45('congctrl', [0 500], [w0 b0], [], N, rho, c);

figure(1); clf; subplot(221);
plot(...
  t1, y1(:,1), 'k', ...
  t1, y1(:,2), 'k', ...
  t1, y1(:,3), 'k',  ...
  t1, y1(:,4), 'k', ...
  t1, y1(:,5), 'k', ...
  t1, y1(:,6), 'k' ...
)
hold on; xy = plot(t1, y1(:,7)/20, 'b');
set(xy, 'LineWidth', AM_data_linewidth);
amaxis();

% Now change the number of sources and rerun from the old initial condition
N = 40;
[weq, beq] = congctrl_equil(N, rho, c)

T = length(t1);
x0 = y1(T, 1:N/10);
b0 = y1(T, 7);
[t2, y2] = ode45('congctrl', [500 1000], [x0 b0], [], N, rho, c);

hold on;
plot(...
  t2, y2(:,1), 'k', ...
  t2, y2(:,2), 'k', ...
  t2, y2(:,3), 'k', ...
  t2, y2(:,4), 'k' ...
)
hold on; xy = plot(t2, y2(:,5)/20, 'b');
set(xy, 'LineWidth', AM_data_linewidth);

xlabel('Time {\itt} [ms]');
ylabel('States {\itw}_{\iti} [pkts/ms], {\itb} [pkts]', 'Rotation', 90);

amtext(800, 11, '{\itb}');
amtext(250, 3, '{\itw}_1-{\itw}_{60}');
amtext(750, 4, '{\itw}_1-{\itw}_{40}');

amprint('congctrl-tcpsim.eps');
