#
# Tumor cell density in glioblastomas
#
# Delete previous workspaces
  rm(list=ls(all=TRUE))
#
# Access ODE integrator
  library("deSolve");
#
# Access functions for numerical solutions
  setwd("f:/comp3/wiley_download/chap9");
  source("pde_1.R");
  source("dss004.R");
  source("dss044.R");
#
# Level of output
#
#   ip = 1 - graphical (plotted) solutions 
#            (u(r,t)) only
#
#   ip = 2 - numerical and graphical output
#
  ip=2;
#
# Select case
  ncase=1;
#
# Parameters
  rho=0.012;umax=62.5;eps=0.01;k=1;d=1;
#
# Grid in r
  nr=101;r0=25;
  r=seq(from=0,to=r0,by=(r0-0)/(nr-1));
#
# Diffusivity
  D=rep(0,nr);
  for(i in 1:nr){
    if(i<=16){D[i]=0.13;}
    if((i>16)&(i<=86)){D[i]=0.65;}
    if(i> 86){D[i]=0.13;}
  }
#
# Display D(r)
  if(ip==2){
    for(i in 1:nr){
      cat(sprintf("\n i = %3d  r = %4.2f  D(r) = %5.3f",
                   i,r[i],D[i]));
    }
  }
#
# dD/dr
  dDdr=rep(0,nr);
  for(i in 1:nr){
    if(i< 16){dDdr[i]=0;}
    if(i==16){dDdr[i]=(D[i+1]-D[i-1])/0.5;}
    if((i>16)&(i<86)){dDdr[i]=0;}
    if(i==86){dDdr[i]=(D[i+1]-D[i-1])/0.5;}
    if(i> 86){dDdr[i]=0;}
  }
#
# Display dD/dr
  if(ip==2){
    for(i in 1:nr){
      cat(sprintf("\n i = %3d  r = %4.2f  dD/dr = %5.3f",
                   i,r[i],dDdr[i]));
    }
  }
#
# Initial conditions
  u0=rep(0,nr);
  fact=1/((2*pi)^0.5*eps);
  for(i in 1:nr){
      u0[i]=exp(-0.5*(r[i]-12.5)^2/eps);
  }
  u0=fact*u0;
  ncall=0;
#
# Write selected parameters
  cat(sprintf("\n\n ncase = %2d",ncase));
#
# Write heading
  if(ip==1){
    cat(sprintf("\n Graphical output only\n"));
  } 
#
# Independent variable for ODE integration
  nout=6;t0=0;tf=15;
  tout=seq(from=t0,to=tf,by=(tf-t0)/(nout-1));
#
# ODE integration
  out=lsodes(y=u0,times=tout,func=pde_1,parms=NULL)
  nrow(out)
  ncol(out)
#
# Arrays for plotting numerical solution
  u_plot1=matrix(0,nrow=nr,ncol=nout);
  u_plot2=matrix(0,nrow=nr,ncol=nout-1);
  for(it in 1:nout){
    for(i in 1:nr){
       u_plot1[i,it]=out[it,i+1];
       if(it>1){u_plot2[i,it-1]=u_plot1[i,it];}
    }
  }
#
# Display numerical solution
  if(ip==2){
    for(it in 1:nout){
      cat(sprintf(
      "\n    t       r    u(r,t)\n"));
      for(i in 1:nr){
        cat(sprintf("%5.1f%8.2f%10.3f\n",
        tout[it],r[i],u_plot1[i,it]));
      }
    }
  }
#
# Calls to ODE routine
  cat(sprintf("\n\n ncall = %5d\n\n",ncall));
#
# Plot u with t = 0
  par(mfrow=c(1,1));    
  matplot(x=r,y=u_plot1,type="l",xlab="r",
          ylab="u(r,t), t=0,3,...,15",xlim=c(0,r0),lty=1,
          main="u(r,t); t=0,3,...,15;",lwd=2,col="black");
#
# Plot u without t = 0
  par(mfrow=c(1,1));    
  matplot(x=r,y=u_plot2,type="l",xlab="r",
          ylab="u(r,t), t=3,...,15",xlim=c(0,r0),lty=1,
          main="u(r,t); t=3,...,15;",lwd=2,col="black");