% plot simulation result
global Dose r k13 k43 cl dr br k loadmax Ntotal phi rho nvar tmax tend lambda conc
% plot 10 mgm-3 data
% data for comparison TiO2 (10 mgm-3)

t10=[0, 12, 26, 54, 82, 89, 110, 117, 236, 341, 579, 684];
b10data=[0, 0.897, 1.674, 2.672, 3.797, 3.99, 3.0, 4.197, 5.62, 7.52, 18.2, 19.6];
l10data=[0,0,0,0,0,0.018,0.047,0.058,0.169,0.762,3.85,5.21];
b10=[3.99, 3.00, 5.62, 7.52, 18.2, 19.6];
l10=[0.018,0.047,0.169,0.762,3.85,5.21];


tss=[89, 110, 236, 341, 579, 684];
nb=[24, 6, 6, 6, 5, 6];
nl=[12, 4, 6, 6, 5, 6];
std1=[0.78, 0.26, 1.2, 0.77, 7.1, 4.0];
std2=[0.004, 0.013, 0.011, 0.018, 1.1, 1.4];


for i=1:6
b10l(i)=b10(i) - 2*(std1(i)/sqrt(nb(i)));
b10h(i)=b10(i) + 2*(std1(i)/sqrt(nb(i)));
l10l(i)=l10(i) - 2*(std2(i)/sqrt(nl(i)));
l10h(i)=l10(i) + 2*(std2(i)/sqrt(nl(i)));
end;

% TiO2 (30 mgm-3)
  t30=[0, 12, 26, 40, 54, 82, 117];
  b30data=[0, 3.228, 3.091, 8.336, 8.354, 10.06, 18.373];
  l30data=[0, 0, 0.013, 0.025, 0.027, 0.168, 1.16];
% TiO2 (50 mgm-3)
  t50=[0,26,54,110,220];
  b50data=[0,6.23,13.3,33.6,57.9];
  l50data=[0,0.134,0.202,2.16,16.3];
% TiO2 (90 mgm-3)
  t90=[0, 12, 26, 40, 54, 82, 117];
  b90data=[0, 11.945, 13.094, 24.00, 31.267, 39.6, 55.634];
  l90data=[0, 0.06, 0.11, 0.360, 0.828, 2.209, 3.17];

% short inhalation
  tt1=[10, 38]; 
  pb1=[0.910,0.476];
  tt2=[10, 38];
  pb2=[1.684,0.998];
  tt3=[10,38];
  pb3=[1.649,1.661];
  tt4=[10, 38];
  pb4=[3.328,2.970];
  tt5=[10,38];
  pb5=[3.875,2.599];
  tt8=[10,38];
  pb8=[3.875,2.599];
  pl8=[21, 33];

% extended inhalation 1
  tt6=[38,185];
  pb6=[2.780,1.200];
  nrb6=[24, 32];
  nrl6=[12, 32];
  stdb6=[0.82, 0.4];
  stdl6=[0.014, 0.071];
  pl6=[0.029,0.119];
  for i=1:2
  pb6l(i)=pb6(i) - 2*(stdb6(i)/sqrt(nrb6(i)));
  pb6h(i)=pb6(i) + 2*(stdb6(i)/sqrt(nrb6(i)));
  pl6l(i)=pl6(i) - 2*(stdl6(i)/sqrt(nrl6(i)));
  pl6h(i)=pl6(i) + 2*(stdl6(i)/sqrt(nrl6(i)));
  end;

% 2
  tt7=[115,227,465];
  pb7=[1.020,0.773,0.264];
  pl7=[0.094, 0.212,0.063];
  nrb7=[6, 6, 5];
  nrl7=[6, 6, 6];
  stdb7=[0.13, 0.14, 0.09];
  stdl7=[0.011, 0.044, 0.058];
  for i=1:3
  pb7l(i)=pb7(i) - 2*(stdb7(i)/sqrt(nrb7(i)));
  pb7h(i)=pb7(i) + 2*(stdb7(i)/sqrt(nrb7(i)));
  pl7l(i)=pl7(i) - 2*(stdl7(i)/sqrt(nrl7(i)));
  pl7h(i)=pl7(i) + 2*(stdl7(i)/sqrt(nrl7(i)));
  end;
  
% 5
  pb10=[13.7];
  stdb10=[2.4];
  nrb10=[6];
  pb10l(1)=pb10(1) - 2*(stdb10(1)/sqrt(nrb10(1)));
  pb10h(1)=pb10(1) + 2*(stdb10(1)/sqrt(nrb10(1)));
  tt10=[42];
% no data for lymph


% 6
  tt11=[38,73];
  pb11=[19.5,9.65];
  pl11=[3.6,5.67];
  stdb11=[4.8, 1.7];
  stdl11=[0.39, 2.7];
  nrb11=[6, 5];
  nrl11=[4, 4];
  for i=1:2
  pb11l(i)=pb11(i) - 2*(stdb11(i)/sqrt(nrb11(i)));
  pb11h(i)=pb11(i) + 2*(stdb11(i)/sqrt(nrb11(i)));
  pl11l(i)=pl11(i) - 2*(stdl11(i)/sqrt(nrl11(i)));
  pl11h(i)=pl11(i) + 2*(stdl11(i)/sqrt(nrl11(i)));
  end;

 % 3 
  tt12=[42];
  nrb12=[5];
  nrl12=[6];
  pb12=[4.790];
  pl12=[3.340];
  stdb12=[1.6];
  stdl12=[1.2];
  pb12l(1)=pb12(1) - 2*(stdb12(1)/sqrt(nrb12(1)));
  pb12h(1)=pb12(1) + 2*(stdb12(1)/sqrt(nrb12(1)));
  pl12l(1)=pl12(1) - 2*(stdl12(1)/sqrt(nrl12(1)));
  pl12h(1)=pl12(1) + 2*(stdl12(1)/sqrt(nrl12(1)));

% 4
  tt13=[38];
  pb13=[6.13];
  pl13=[1.07];
  nrb13=[6];
  nrl13=[6];
  stdb13=[1.6];
  stdl13=[0.71];
  pb13l(1)=pb13(1) - 2*(stdb13(1)/sqrt(nrb13(1)));
  pb13h(1)=pb13(1) + 2*(stdb13(1)/sqrt(nrb13(1)));
  pl13l(1)=pl13(1) - 2*(stdl13(1)/sqrt(nrl13(1)));
  pl13h(1)=pl13(1) + 2*(stdl13(1)/sqrt(nrl13(1)));

  alvburden=(100/123)*(b10data);
  instburden=(23/100)*alvburden;


figure(1)

subplot(221),plot(time10,burden10,t10,b10data,'ko',tss,b10l,'.',tss,b10h,'.')
axis([0 690 0 30])
title(' TiO_2 - 10 mg/m^3')
xlabel('Time - (days)')
ylabel('Lung burden (mg)')

subplot(222),plot(time10,lburden10,t10,l10data,'ko',tss,l10l,'.',tss,l10h,'.')
axis([0 690 0 10.0])
xlabel('Time - (days)')
ylabel('Lymph burden (mg)')

subplot(223), plot(time10, PMN)
axis([0 690 0 4])
xlabel('Time - (days)')
ylabel(' N^o of PMN (x10^6)')

subplot(224), plot(burden10, PMN)
axis([0 20 0 4])
xlabel('Lung burden (mg)')
ylabel('N^o of PMN (x10^6)')