%% .. Polymerization
% ...  
% ... 13.03.2014:
% ...
% ... Milan Bernauer, ICT, Prague 2014
%
clear all; clc;
global kd kp kt n nm
kd = 1.45e-5; 
kp = 4.4e2;
kt = 1.2e6;
cI = 1.0e-2;
cM = 3.0;

% --- Set cP(n)
nm = [50 150 500 1000];
n  = size(nm,2);

% --- Set initial conditions
c0(1) = cI;
c0(2) = cM;
c0(3:n+2) = 0.0;

% --- Set time interval [s]
tr = 10*3.6e3;
% --- Set steps (to be compatible with OCTAVE ode15s wrapper)
incr = tr/50;

% --- Calling ODE solver
[t,c] = ode15s(@dcdt,[0:incr:tr],c0);

% --- Plot results (two y axes)
[AX,H1,H2] = plotyy(t/3600,c(:,1:2), ...
    t/3600,c(:,3:n+2));

title('Polymererization - Steady state hypotesis')
grid on
xlabel('Time / h')
ylabel('c / (mol l^{-1})')
set(get(AX(1),'Ylabel'),'String','c_I, c_M') 
set(get(AX(2),'Ylabel'),'String','c_P(n)')
set(H1,'LineStyle','--','LineWidth',2)
set(H2,'LineStyle','-','LineWidth',2)

% --- Set labels (legend)
leg = num2str(nm'); % Will not work in OCTAVE

legend(H1,'c_I','c_M','Location','northwest');
legend(H2,leg,'Location','north');
