function Figure6_7

% Matlab code for Figure 6.7 in book:
% Radek Erban and Jon Chapman, Stochastic Modelling of Reaction-Diffusion
% Processes, Cambridge University Press, 2019

close all;


D=0.0001;
h=0.01;
xdet=[h/2:h:1-h/2];
w=size(xdet);
boxdeterm=w(2);
dt=0.01;
Diff=D/(h*h);
detden=zeros(boxdeterm,1);
influxrate=0.3125;
k1=0.001;

for t=0:dt:30*60
    detdenold=detden;
    ii=1;     
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)-detdenold(ii))-dt*k1*detdenold(ii);
    for ii=2:boxdeterm-1
       detden(ii)=detdenold(ii)+dt*Diff*(detdenold(ii+1)+detdenold(ii-1)-2*detdenold(ii))-dt*k1*detdenold(ii);
    end  
    ii=boxdeterm;
       detden(ii)=detdenold(ii)+dt*Diff*(2*detdenold(ii-1)-2*detdenold(ii))-dt*k1*detdenold(ii);

    for ii=1:boxdeterm/5
       detden(ii)=detden(ii)+dt*influxrate;
    end  
end

bw=26;
wr=54;

detblueX=[0; 0; xdet(1:bw)'; xdet(bw); 0];
detblueY=[0; detden(1); detden(1:bw); 0; 0];

detwhiteX=[xdet(bw); xdet(bw:wr)'; xdet(wr); xdet(bw)];
detwhiteY=[0; detden(bw:wr); 0; 0];

detredX=[xdet(wr); xdet(wr:boxdeterm)'; 1; 1];
detredY=[0; detden(wr:boxdeterm); detden(boxdeterm); 0];

figure(1);
set(gca,'Fontsize',18);
rectangle('Position',[0,0,1,155],'FaceColor',[0.86 0.86 0.86]);
hold on;
hh1=bar([2],[0],'b');
hh2=bar([2],[0],'w');
hh3=bar([2],[0],'r');
plot([0 xdet 1],[detden(1); detden; detden(boxdeterm)],'k','Linewidth',4);
fill(detblueX,detblueY,'b');
fill(detwhiteX,detwhiteY,'w');
fill(detredX,detredY,'r');
h=line([xdet(wr) 1],[detden(wr) detden(wr)]);
set(h,'Color','k','LineStyle','--');
h=line([xdet(bw) 1],[detden(bw) detden(bw)]);
set(h,'Color','k','LineStyle','--');
xlabel('$x$ [mm]','interpreter','latex');
ylabel('number of molecules','interpreter','latex');
axis([0 1 0 155]);
set(gca,'YTick',[0 20 40 60 80 100 120 140]);
hh=legend([hh1,hh2,hh3],'$N > 80$','$80 \geq N > 30$','$30 \geq N$');
set(hh,'Fontsize',18,'interpreter','latex');
box on;
set(gca,'Fontsize',18);

histogram=load('data_Figure6_1.dat');
ww=size(histogram);
boxes=ww(1);
mesh=[1/(2*boxes):1/boxes:1-1/(2*boxes)]';

morblue=zeros(boxes,1);
morwhite=zeros(boxes,1);
morred=zeros(boxes,1);

for ii=1:boxes
    if (histogram(ii,2)>80)
       morblue(ii)=histogram(ii,2);
    else
       if (histogram(ii,2)>30)
           morwhite(ii)=histogram(ii,2);
       else
          morred(ii)=histogram(ii,2);
       end
    end
end

figure(2);
set(gca,'Fontsize',18);
rectangle('Position',[0,0,1,155],'FaceColor',[0.86 0.86 0.86]);
hold on;
hh1=bar(mesh,[morblue],'b');
hh2=bar(mesh,[morwhite],'w');
hh3=bar(mesh,[morred],'r');
h=line([0 1],[30 30]);
set(h,'Color','k','LineStyle','--');
h=line([0 1],[80 80]);
set(h,'Color','k','LineStyle','--');
xlabel('$x$ [mm]','interpreter','latex');
ylabel('number of molecules','interpreter','latex');
axis([0 1 0 155]);
set(gca,'YTick',[0 20 40 60 80 100 120 140]);
hh=legend([hh1,hh2,hh3],'$N > 80$','$80 \geq N > 30$','$30 \geq N$');
set(hh,'Fontsize',18,'interpreter','latex');
box on;
set(gca,'Fontsize',18);

% ======================================
% Turing Patterns
% ======================================

boxes=40;

load data_Figure6_7_c_d.dat;
dataSTO=data_Figure6_7_c_d;

morgreen=zeros(boxes,1);
morblack=zeros(boxes,1);
morgreen2=zeros(boxes,1);
morblack2=zeros(boxes,1);

for ii=1:boxes
    if (dataSTO(ii,2)>200)
       morgreen(ii)=dataSTO(ii,2);
    else
       morblack(ii)=dataSTO(ii,2);
    end
    if (dataSTO(ii,3)>75)
       morblack2(ii)=dataSTO(ii,3);
    else
       morgreen2(ii)=dataSTO(ii,3);
    end
end

figure(3);
set(gca,'Fontsize',18);
bar(dataSTO(:,1),morgreen(:,1),'g');
hold on;
bar(dataSTO(:,1),morblack(:,1),'k');
xlabel('$x$ [mm]','interpreter','latex');
ylabel('number of $A$ molecules','interpreter','latex');
text(0.027,476,'time=30 min','Fontsize',18,'interpreter','latex');
h=line([0 1],[200 200]);
set(h,'Color','k','LineStyle','--');
axis([0 1 0 510]);
set(gca,'Fontsize',18);

figure(4);
set(gca,'Fontsize',18);
bar(dataSTO(:,1),morgreen2(:,1),'g');
hold on;
bar(dataSTO(:,1),morblack2(:,1),'k');
xlabel('$x$ [mm]','interpreter','latex');
ylabel('number of $B$ molecules','interpreter','latex');
text(0.027,112,'time=30 min','Fontsize',18,'interpreter','latex');
h=line([0 1],[75 75]);
set(h,'Color','k','LineStyle','--');
axis([0 1 0 120]);
set(gca,'Fontsize',18);

