close all; clear; clc;

%% A+B->A'+B' but B'->B as a function of kBp
%
figure(1); clf;
a = subplot(1,2,1); cla; hold on;
a.XScale = 'log'; axis square;
xlim([1e-2 1e2]); ylim([0 1]);
legend off; grid on;
alpha = 0.4; % Transparency of the background curves
clear h;

kBs = 10.^(2:-1:-4);
kBs = 10.^linspace(2,-4,10);

cm = plasma(numel(kBs)+1);
for kBi = 1:numel(kBs)
    y0 = [1 1 0 0]';
    kAB = 1; kBp = kBs(kBi);
    f = @(t,y) kAB*[-1 -1 1 1]'*y(1)*y(2) + [0 1 0 -1]'*y(4)*kBp;
    [t,y] = ode15s(f,[0 1e2],y0);
    
    color = cm(kBi,:);
    h(kBi) = plot(t,y(:,1),'-','color',color,'linewidth',1);
    color = alpha*cm(kBi,:)+(1-alpha)*[1 1 1];
    l = plot(t,y(:,4),'-','color',color,'linewidth',1);
    uistack(l,'bottom');
    drawnow;
    
end
y0 = [1 0];
f = @(t,y) kAB*[-1 0;1 0]*y;
[t,y2] = ode15s(f,t,y0);

plot(t,y2(:,1),'k--','LineWidth',1)
temp = 1./(1+kAB.*t);
plot(t,temp,'k--','LineWidth',1)
%plot(t,1-temp,'k-')
legend(h([1,4,7,end]),'2','0','-2','-4','location','northwest');
grid on; box on;
xlabel('time / s'); ylabel('conc.');

%% A->B B->C+D C->A and D->A
ratios = 10.^linspace(-1,1,7);
y0 = [1 0 0 0]';
kAB = 1000; kBCD = 1000; kCA = 1;
figure(1);% clf;
subplot(1,2,2); cla; axis square; hold on;
legend('off');
cm = lines;
for j = 1:numel(ratios)
    kDA = kCA*ratios(j);
    K = [-kAB   0      +kCA +kDA ;
        +kAB  -kBCD    0    0   ;
        0    +kBCD/2 -kCA  0   ;
        0    +kBCD/2  0   -kDA];
    f = @(t,y) K*y;
    [t,y] = ode15s(f,[0 1e2],y0);
    
    style = '-';
    alpha = ((j-1)/numel(ratios)); alpha = alpha.^(0.8);%alpha = sqrt(alpha);
    color = cm(1,:).*alpha + (1-alpha).*[1 1 1];
    color = color.*(0.4) + [1 1 1]*0.6;
    h1 = plot(t,y(:,1),'color',color,'linewidth',1,'linestyle',style);
    color = cm(2,:).*alpha + (1-alpha).*[1 1 1];
    color = color.*(0.4) + [1 1 1]*0.6;
    h2 = plot(t,y(:,2),'color',color,'linewidth',1,'linestyle',style);
    color = cm(3,:).*alpha + (1-alpha).*[1 1 1];
    h3 = plot(t,y(:,3),'color',color,'linewidth',1,'linestyle',style);
    color = cm(4,:).*alpha + (1-alpha).*[1 1 1];
    h4 = plot(t,y(:,4),'color',color,'linewidth',1,'linestyle',style);
    if j == numel(ratios)
        h = [h1 h2 h3 h4];
    end
    a = gca; a.XScale = 'log';
    grid on; box on;
    xlim([1e-5 1e2]);
    drawnow;
end

kDA = kCA;
K = [-kAB   0      +kCA +kDA ;
    +kAB  -kBCD    0    0   ;
    0    +kBCD/2 -kCA  0   ;
    0    +kBCD/2  0   -kDA];
f = @(t,y) K*y;
[t,y] = ode15s(f,[0 1e2],y0);
h(5) = plot(t,y(:,3),'k--','linewidth',1);
legend(h,'[A]','[B]','[C]','[D]','r=1','location','northwest');
ylabel('conc.');
xlabel('time / s');
