%
% Four-section retinal O2 transport model equation
%
% Clear previous files
  clear all
  clc
%
% Model parameters shared with other routines
%
% Model parameters
  global   nir     nor     nfl     ncc...
         zl_ir   zl_or   zl_fl   zl_cc...
         zg_ir   zg_or   zg_fl   zg_cc...
          D1ir    D1or    D1fl    D1cc...
          D2ir    D2or    D2fl    D2cc...          
               k1ir_or k1or_fl k1fl_cc...
               k2ir_or k2or_fl k2fl_cc...               
          k1ir    k1or    k1fl    k1cc...
          k2ir    k2or    k2fl    k2cc...
         u1irt   u1ort   u1flt   u1cct... 
         pir_s   pcc_s   ncall   ncase      nt
%
% Select case
%
%   ncase = 1: VEGF rates = 0
%
%   ncase = 2: VEGF rates ne 0
  ncase=2;       
%
% Number of grid points
  nir=11; nor=11; nfl=11; ncc=11;
  nt=nir+nor+nfl+ncc;
%
% Lengths of four sections (microns)
  zl_ir=200; zl_or=200; zl_fl=200; zl_cc=200;
%
% Spatial grids;
  zg_ir=[0:zl_ir/10:zl_ir]'; zg_or=[0:zl_or/10:zl_or]';
  zg_fl=[0:zl_fl/10:zl_fl]'; zg_cc=[0:zl_cc/10:zl_cc]';
%
% Diffusivities (microns^2/sec)
  D1ir=1.0e+04; D1or=1.0e+04; D1fl=1.0e+04; D1cc=1.0e+04;
  D2ir=1.0e+03; D2or=1.0e+03; D2fl=1.0e+03; D2cc=1.0e+03;
%
% Interface coefficients
  k1ir_or=1; k1or_fl=1; k1fl_cc=1;
  k2ir_or=1; k2or_fl=1; k2fl_cc=1;  
%
% Metabolism rates
  k1ir=0.1; k1or=0.1; k1fl=0.1; k1cc=0.1;
%
% VEGF rates
  if(ncase==1) k2ir=0      ; k2or=0      ; k2fl=0      ; k2cc=0      ; end 
  if(ncase==2) k2ir=1.0e-05; k2or=1.0e-05; k2fl=1.0e-05; k2cc=1.0e-05; end   
%
% Normal breathing O2 concentration (mm hg)
  pir_s= 20;
  pcc_s=100;  
  pcc_s=  0;
%
% O2 thresholds
  u1irt=10; u1ort=10; u1flt=10; u1cct=10;
%
% Independent variable t
  t0=0.0; tf=7200; tout=[t0:tf/6:tf]'; 
  nout=7;
%  
% 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;
    options=odeset(options,'JPattern',S)
    [t,u]=ode15s(@pde_1,tout,u0,options);   
  end 
%
% One vector to nine vectors
%
% u1
    for it=1:nout
      for i=1:nir u1ir(it,i)=u(it,i);                end
      for i=1:nor u1or(it,i)=u(it,i+nir);            end
      for i=1:nfl u1fl(it,i)=u(it,i+nir+nor);        end
      for i=1:ncc u1cc(it,i)=u(it,i+nir+nor+nfl);    end
      if(it>=2)
        u1ir(it,1)  =pir_s;           
        u1cc(it,ncc)=pcc_s; 
      end
    end
%
% u2
    for it=1:nout
      for i=1:nir u2ir(it,i)=u(it,nt+i);             end
      for i=1:nor u2or(it,i)=u(it,nt+i+nir);         end
      for i=1:nfl u2fl(it,i)=u(it,nt+i+nir+nor);     end
      for i=1:ncc u2cc(it,i)=u(it,nt+i+nir+nor+nfl); end
    end
%
%   Display tabulated numerical solution
%
%   Inner retina
    for it=1:nout
      fprintf('\n    t     zg_ir   u1ir(z,t)')
      fprintf('\n    t     zg_ir   u2ir(z,t)\n')
      for i=1:nir
        fprintf('%5.0f%10.0f%12.2f\n'  ,t(it)/60,zg_ir(i),u1ir(it,i))
        fprintf('%5.0f%10.0f%12.2f\n\n',t(it)/60,zg_ir(i),u2ir(it,i))        
      end 
%
%     Outer retina
      fprintf('\n    t     zg_or   u1or(z,t)')
      fprintf('\n    t     zg_or   u2or(z,t)\n')
      for i=1:nor
        fprintf('%5.0f%10.0f%12.2f\n'  ,t(it)/60,zg_or(i),u1or(it,i))
        fprintf('%5.0f%10.0f%12.2f\n\n',t(it)/60,zg_or(i),u2or(it,i))
      end  
%
%     Fluid layer      
      fprintf('\n    t     zg_fl   u1fl(z,t)\n')
      fprintf('\n    t     zg_or   u2fl(z,t)\n')      
      for i=1:nfl
        fprintf('%5.0f%10.0f%12.2f\n'  ,t(it)/60,zg_fl(i),u1fl(it,i))
        fprintf('%5.0f%10.2f%12.2f\n\n',t(it)/60,zg_fl(i),u2fl(it,i))        
      end  
%
%     Choriocapillaris      
      fprintf('\n    t     zg_cc   u1cc(z,t)\n')
      for i=1:ncc
        fprintf('%5.0f%10.0f%12.2f\n'  ,t(it)/60,zg_cc(i),u1cc(it,i))
        fprintf('%5.0f%10.2f%12.2f\n\n',t(it)/60,zg_cc(i),u2cc(it,i))        
      end  
%
%   Next t
    end   
    fprintf('\n    ncall = %4d\n',ncall);
%
%   Four plots for u1ir, u1or, u1fl, u1cc
    figure(2)
    subplot(2,2,1)
    plot(zg_ir,u1ir,'-')
    axis([0 200 0 100])
    xlabel('zg-ir, \mum'); ylabel('u1ir(z,t), mm hg O_2');
    title('Inner retina O_2, t=0,20,...,120 min');  
    subplot(2,2,2)
    plot(zg_or,u1or,'-')
    axis([0 200 0 100])
    xlabel('zg-or, \mum'); ylabel('u1or(z,t), mm hg O_2'); 
    title('Outer retina O_2, t=0,20,...,120 min');  
    subplot(2,2,3)
    plot(zg_fl,u1fl,'-')
    axis([0 200 0 100])
    xlabel('zg-fl, \mum'); ylabel('u1fl(z,t), mm hg O_2');
    title('Fluid layer O_2, t=0,20,...,120 min');
    subplot(2,2,4)
    plot(zg_cc,u1cc,'-')
    axis([0 200 0 100])
    xlabel('zg-cc, \mum'); ylabel('u1cc(z,t), mm hg O_2');
    title('Choriocapillaris O_2, t=0,20,...,120 min');  
%
%   Four plots for u2ir, u2or, u2fl, u2cc
    figure(3)
    subplot(2,2,1)
    plot(zg_ir,u2ir,'-')
    axis([0 200 0 1])
    xlabel('zg-ir, \mum'); ylabel('u2ir(z,t), ng/ml');
    title('Inner retina VEGF, t=0,20,...,120 min');  
    subplot(2,2,2)
    plot(zg_or,u2or,'-')
    axis([0 200 0 1])
    xlabel('zg-or, \mum'); ylabel('u2or(z,t), ng/ml'); 
    title('Outer retina VEGF, t=0,20,...,120 min');  
    subplot(2,2,3)
    plot(zg_fl,u2fl,'-')
    axis([0 200 0 1])
    xlabel('zg-fl, \mum'); ylabel('u2fl(z,t), ng/ml');
    title('Fluid layer VEGF, t=0,20,...,120 min');
    subplot(2,2,4)
    plot(zg_cc,u2cc,'-')
    axis([0 200 0 1])
    xlabel('zg-cc, \mum'); ylabel('u2cc(z,t), ng/ml');
    title('Choriocapillaris VEGF, t=0,20,...,120 min');  
   
    
    
    

    
    
