%
% Three-PDE ATG model
%
% Clear previous files
  clear all
  clc
%
% Model parameters shared with other routines
  global    Nn      Nt      Ch...
            rl      ru       r       n...
           rn1     rn2      Kn...
           rt1      Dt      Kt...
           rh1     rh2      Dh...
         ncall   ncase  
%
% Case
%
%   ncase = 1 - four 2D plots of dependent variables Nn, Nt, Ch, pH vs t
%
%   ncase = 2 - four 3D plots of dependent variables Nn, Nt, Ch, pH vs t 
%
%   ncase = 3 - case = 1 with 17 supplemental plots for the derivatives 
%               in t and the individual PDE RHS terms
  ncase=3;
%
% Radial grid
  rl=0; ru=0.5; n=101;
  for i=1:n
    r(i)=rl+(ru-rl)*(i-1)/(n-1);
  end  
% 
% Nn PDE
  rn1=1.0e-06; rn2=1; Kn=5.0e+07;
% 
% Nt PDE
  rt1=1.0e-06; Dt=2.0e-10; Kt=5.0e+07;
% 
% Ch PDE
  rh1=2.2e-17; rh2=1.1e-04; Dh=5.0e-06;  
%
% Independent variable t for ncase = 1, 3
  if(ncase==1 | ncase==3)
    t0=0.0;
    tf=5.0e+06;
    tout=[t0:1.0e+06:tf]';
    nout=6;
  end 
%
% Independent variable t for ncase = 2 (to give additional resolution in t)
  if(ncase==2)
    t0=0.0;
    tf=5.0e+06;
    tout=[t0:0.25e+06:tf]';
    nout=21;
  end    
%  
% Initial condition 
  u0=inital_1(t0);
%
% ODE integration
  mf=2;
  reltol=1.0e-06; abstol=1.0e-05;
  options=odeset('RelTol',reltol,'AbsTol',abstol);
%
% Explicit (nonstiff) integration
  if(mf==1)[t,u]=ode45(@pde_1,tout,u0,options); end
%
% Implicit (sparse stiff) integration
  if(mf==2)
    S=jpattern_num_1;
%   pause
    options=odeset(options,'JPattern',S)
    [t,u]=ode15s(@pde_1,tout,u0,options);   
  end 
%
% One vector to three vectors
  for it=1:nout
  for i=1:n
    Nn(it,i)=u(it,i);
    Nt(it,i)=u(it,i+n);
    Ch(it,i)=u(it,i+2*n);
  end
  end
%
%   Display selected output
    for it=1:nout
      fprintf('\n     t       r       Nn(it,i)       Nt(it,i)       Ch(it,i)\n');
      for i=1:10:n  
        fprintf('%6.1f%8.3f%15.4e%15.4e%15.4e\n',...
              t(it)/(60*60*24),r(i),Nn(it,i),Nt(it,i),Ch(it,i));
      end
    end  
    fprintf('\n    ncall = %4d\n',ncall);
%
% 2D plots
%
% ncase = 1, 3: 2D Nn, Nt, Ch, pH plots in t (figures 1 - 4)
  if(ncase==1 | ncase==3)
    pH=-log10(4.0e-08+Ch);
    figure(1)
    plot(r,Nn,'-'); xlabel('r (cm)'); ylabel('Nn (cells/cm^3)')
    title('Nn (normal); 3-PDE ATG model')
    figure(2)
    plot(r,Nt,'-'); xlabel('r (cm)'); ylabel('Nt (cells/cm^3)')
    title('Nt (tumor); 3-PDE ATG model')
    figure(3)
    plot(r,Ch,'-'); xlabel('r (cm)'); ylabel('Ch (M = gm mols H^+/liter)')
    title('Ch (excess H^+); 3-PDE ATG model')
    figure(4)
    plot(r,pH,'-'); xlabel('r (cm)'); ylabel('pH')
    title('pH; 3-PDE ATG model')
  end  
%
% ncase = 2: 3D Nn, Nt, Ch, pH plots in t (figures 1 - 4)
  if(ncase==2)
    tp=t/(60*60*24);  
    pH=-log10(4.0e-08+Ch);
    figure(1)
    surf(r,tp,Nn)
    view(-40,50);
    xlabel('r (cm)'); ylabel('t (days)'); 
    zlabel('Nt (normal cells/cm^3)');
    figure(2)
    surf(r,tp,Nt)
    view(-25,50)
    xlabel('r (cm)'); ylabel('t (days)'); 
    zlabel('Nt (tumor cells/cm^3)');
    figure(3)
    surf(r,tp,Ch)
    xlabel('r (cm)'); ylabel('t (days)'); 
    zlabel('Ch (excess gm mols H^+/liter)');
    figure(4)
    surf(r,tp,pH)
    view(-37.5,60);
    xlabel('r (cm)'); ylabel('t (days)'); 
    zlabel('pH');
  end  
%
% case = 3 supplemental 2D plots for derivatives in t (figures 5 - 7)
  if(ncase==3)
    for it=1:nout
      u1d=u(it,:);  
      ut=pde_1(t(it),u1d);
      Nnt(it,:)=ut(1:n);
      Ntt(it,:)=ut(n+1:2*n);
      Cht(it,:)=ut(2*n+1:3*n);
    end  
    figure(5)
    plot(r,Nnt,'-'); xlabel('r (cm)'); ylabel('Nnt (cells/cm^3-s)')
    title('Nnt (normal); 3-PDE ATG model') 
    figure(6)
    plot(r,Ntt,'-'); xlabel('r (cm)'); ylabel('Ntt (cells/cm^3-s)')
    title('Ntt (normal); 3-PDE ATG model')
    figure(7)
    plot(r,Cht,'-'); xlabel('r (cm)'); ylabel('Cht (M/s = gm mols H^+/liter-s)')
    title('Cht (excess H^+); 3-PDE ATG model') 
%
% Supplemental 2D plots for RHS terms in derivatives in t (figures 8 - 21)
    for it=1:nout
      u1d=u(it,:);  
      [ut,term]=pde_1c(t(it),u1d);
      RHS1(it,:)=term(1,1:n);
      RHS2(it,:)=term(2,1:n);
      RHS3(it,:)=term(3,1:n);
      RHS4(it,:)=term(4,1:n);
      RHS5(it,:)=term(5,1:n);
      RHS6(it,:)=term(6,1:n);
      RHS7(it,:)=term(7,1:n);
      RHS8(it,:)=term(8,1:n);
    end 
%
%   Nn PDE
    figure(8)
    plot(r,RHS1,'-'); xlabel('r (cm)'); ylabel('Nn PDE RHS term 1')
    title('Nn PDE RHS term 1; Nn PDE') 
    figure(9)
    plot(r,RHS2,'-'); xlabel('r (cm)'); ylabel('Nn PDE RHS term 2')
    title('Nn PDE RHS term 2; Nn PDE') 
    figure(10)
    plot(r,RHS1(1,:),'+',r,RHS2(1,:),'o',r,Nnt(1,:),'-'); 
    xlabel('r (cm)'); ylabel('Nn PDE RHS terms 1,2,sum; t=t0')
    title('term 1 - +, term 2 - o, sum = solid; t = t0; Nn PDE')
    figure(11)
    plot(r,RHS1(nout,:),'+',r,RHS2(nout,:),'o',r,Nnt(nout,:),'-'); 
    xlabel('r (cm)'); ylabel('Nn PDE RHS terms 1,2,sum; t=tf')
    title('term 1 - +, term 2 - o, sum = solid; t = tf; Nn PDE')
%
%   Nt PDE
    figure(12)
    plot(r,RHS3,'-'); xlabel('r (cm)'); ylabel('Nt PDE RHS term 1')
    title('Nt PDE RHS term 1; Nt PDE') 
    figure(13)
    plot(r,RHS4,'-'); xlabel('r (cm)'); ylabel('Nt PDE RHS term 2')
    title('Nt PDE RHS term 2; Nt PDE') 
    figure(14)
    plot(r,RHS5,'-'); xlabel('r (cm)'); ylabel('Nt PDE RHS term 3')
    title('Nt PDE RHS term 3; Nt PDE') 
    figure(15)
    plot(r,RHS3(1,:),'+',r,RHS4(1,:),'o',r,RHS5(1,:),'x',r,Ntt(1,:),'-'); 
    xlabel('r (cm)'); ylabel('Nt PDE RHS terms 1,2,3,sum; t = t0')
    title('term 1 - +, term 2 - o, term 3 = x; sum = solid; t = t0; Nt PDE')
    figure(16)
    plot(r,RHS3(nout,:),'+',r,RHS4(nout,:),'o',r,RHS5(nout,:),'x',r,Ntt(nout,:),'-'); 
    xlabel('r (cm)'); ylabel('Nt PDE RHS terms 1,2,sum; t = tf')
    title('term 1 - +, term 2 - o, term 3 - x, sum = solid; t = tf; Nt PDE')
%
%   Ch PDE
    figure(17)
    plot(r,RHS6,'-'); xlabel('r (cm)'); ylabel('Ch PDE RHS term 1')
    title('Ch PDE RHS term 1; Ch PDE') 
    figure(18)
    plot(r,RHS7,'-'); xlabel('r (cm)'); ylabel('Ch PDE RHS term 2')
    title('Ch PDE RHS term 2; Ch PDE') 
    figure(19)
    plot(r,RHS8,'-'); xlabel('r (cm)'); ylabel('Ch PDE RHS term 3')
    title('Ch PDE RHS term 3; Ch PDE') 
    figure(20)
    plot(r,RHS6(1,:),'+',r,RHS7(1,:),'o',r,RHS8(1,:),'x',r,Cht(1,:),'-'); 
    xlabel('r (cm)'); ylabel('Ch PDE RHS terms 1,2,3,sum; t = t0')
    title('term 1 - +, term 2 - o, term 3 = x; sum = solid; t = t0; Ch PDE')
    figure(21)
    plot(r,RHS6(nout,:),'+',r,RHS7(nout,:),'o',r,RHS8(nout,:),'x',r,Cht(nout,:),'-'); 
    xlabel('r (cm)'); ylabel('Ch PDE RHS terms 1,2,sum; t = tf')
    title('term 1 - +, term 2 - o, term 3 - x, sum = solid; t = tf; Ch PDE')
  end    
    

    


