#
# 2-ODE patient, 8-PDE artificial liver support 
# system model
#
# Delete previous workspaces
  rm(list=ls(all=TRUE))
#
# Access ODE integrator
  library("deSolve");
#
# Access functions for numerical solution
  setwd("f:/comp3/wiley_download/chap6");
  source("pde_1.R");
  source("dss020.R");
  source("dss044.R");
#
# Level of output
#
#   ip = 1 - graphical (plotted) solutions 
#            (u1(x,t), u2(x,t)) only
#
#   ip = 2 - numerical and graphical solutions
#
  ip=2;
#
# Parameters
#
# 2-ode model
    VP=5000; VE=(5/3)*VP; lam12=0.0047; lam21=0.0017; QB=150;
   QE=0.005;     alpha=1;    y2in=0.01;        y10=1;  y20=1;
  lam21=0.017;
  cat(sprintf("\n  VP = %4.1f  VE = %4.1f\n",VP,VE));  
  cat(sprintf("\n  lam12 = %8.5f  lam21 = %8.5f\n",lam12,lam21));
  cat(sprintf("\n  QB = %4.1f  QE = %8.5f  alpha = %5.2f\n",
              QB,QE,alpha));
  cat(sprintf("\n  y2in = %5.3f  y10 = %5.3f  y20 = %5.3f\n",
              y2in,y10,y20));
#
# 8-pde model
#
# Initial concentrations
  u10=0;u20=0;u30=0;u40=0;
  u50=0;u60=0;u70=0;u80=0;
#
# u1,u2
  z12=25;v1=-2.5;v2=2.5;
  R1=1;R2=1;k12=0.05;
#
# u3,u4
  z34=25;eps=0.3;k1=1;k2=10;
  v3=-2.5;D3=0.01;k34=0.05;
#
# u5,u6
  z56=25;
  v5=-2.5;D5=0.01;k56=0.05;
#
# u7,u8
  z78=25;v7=-2.5;v8=2.5;u8e=0;
  R7=1;R8=1;k78=0.05;
#
# Independent variable for ODE integration
  nout=37;t0=0;tf=360;
  tout=seq(from=t0,to=tf,by=(tf-t0)/(nout-1));
# 
# Initial conditions
  nz=21;
  u0=rep(0,2+8*nz);
  for(i in 1:nz){
    u0[i]=u10;
    u0[i+nz]=u20;
    u0[i+2*nz]=u30;
    u0[i+3*nz]=u40;
    u0[i+4*nz]=u50;
    u0[i+5*nz]=u60;
    u0[i+6*nz]=u70;
    u0[i+7*nz]=u80;
  }
  u0[1+8*nz]=y10;
  u0[2+8*nz]=y20;
  ncall=0;
#
# ODE integration
  out=lsodes(func=pde_1,y=u0,times=tout,parms=NULL)
#
# Save and display numerical solution
  u1=rep(0,nout);u2=rep(0,nout);u3=rep(0,nout);u4=rep(0,nout);
  u5=rep(0,nout);u6=rep(0,nout);u7=rep(0,nout);u8=rep(0,nout);
  y1=rep(0,nout);y2=rep(0,nout);
  for(it in 1:nout){
    u1[it]=out[it,2]     ;u2[it]=out[it,2*nz+1];
    u3[it]=out[it,2*nz+2];u5[it]=out[it,4*nz+2];
    u4[it]=out[it,3*nz+2];u6[it]=out[it,5*nz+2];
    u7[it]=out[it,6*nz+2];u8[it]=out[it,8*nz+1];
    y1[it]=out[it,8*nz+2];y2[it]=out[it,8*nz+3];
  }
  if(ip==2){
    for(it in 1:nout){
      if(it==1){
        cat(sprintf("\n      t     y1(t)     y2(t)\n"));
      }
    cat(sprintf("%7.1f%10.4f%10.4f\n",tout[it],y1[it],y2[it]));
    }
  }
#
# Calls to ODE routine
  cat(sprintf("\n\n  ncall = %3d\n\n",ncall));
#
# Plot y1, y2
  par(mfrow=c(1,1))
  plot(tout,y1,
     xlab="t (min)",ylab="y1,y2 (micromol/ml)",
     main="y1(t),y2(t)",type="l",lwd=2,
     xlim=c(0,400),ylim=c(0.4,1.3));
   points(tout,y1, pch="1",lwd=2);
    lines(tout,y2,type="l",lwd=2);
   points(tout,y2, pch="2",lwd=2);
#
# Plot u1(z=0,t), u2(z=z12,t)
  par(mfrow=c(1,1))
  plot(tout,u1,
     xlab="t (min)",ylab="u1,u2 (micromol/ml)",
     main="u1(z=0,t),u2(z=z12,t)",type="l",lwd=2,
     xlim=c(0,400));
   points(tout,u1, pch="1",lwd=2);
    lines(tout,u2,type="l",lwd=2);
   points(tout,u2, pch="2",lwd=2);
#
# Plot u3(z=0,t), u5(z=0,t)
  par(mfrow=c(1,1))
  plot(tout,u3,
     xlab="t (min)",ylab="u3,u5 (micromol/ml)",
     main="u3(z=0,t),u5(z=0,t)",type="l",lwd=2,
     xlim=c(0,400));
   points(tout,u3, pch="3",lwd=2);
    lines(tout,u5,type="l",lwd=2);
   points(tout,u5, pch="5",lwd=2);
#
# Plot u4(z=0,t), u6(z=0,t)
  par(mfrow=c(1,1))
  plot(tout,u4,
     xlab="t (min)",ylab="u4,u6 (micromol/ml)",
     main="u4(z=0,t),u6(z=0,t)",type="l",lwd=2,
     xlim=c(0,400));
   points(tout,u4, pch="4",lwd=2);
    lines(tout,u6,type="l",lwd=2);
   points(tout,u6, pch="6",lwd=2);
#
# Plot u7(z=0,t), u8(z=z78,t)
  par(mfrow=c(1,1))
  plot(tout,u7,
     xlab="t (min)",ylab="u7,u8 (micromol/ml)",
     main="u7(z=0,t),u8(z=z78,t)",type="l",lwd=2,
     xlim=c(0,400));
   points(tout,u7, pch="7",lwd=2);
    lines(tout,u8,type="l",lwd=2);
   points(tout,u8, pch="8",lwd=2);