%
% 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...
               k1ir_or k1or_fl k1fl_cc...
          k1ir    k1or    k1fl    k1cc...
         u1ort    k2or                ...
         pir_s   pcc_s   ncall   ncase      nt
%
% Select case
%
%   ncase = 1: Normal O2 at z=z_cc
%
%   ncase = 2: Reduced O2 at z=z_cc
  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/s)
  D1ir=1.0e+04; D1or=1.0e+04; D1fl=1.0e+04; D1cc=1.0e+04;
%
% Interface coefficients
  kb_ir=1; k1ir_or=1; k1or_fl=1; k1fl_cc=1;
%
% Metabolism rates
  k1ir=0.1; k1or=0.1; k1fl=0.1; k1cc=0.1;
%
% Photoreceptor rate
  k2or=2.0e+01;
%
% Normal breathing O2 concentration (mm hg)
  pir_s= 20;
  pcc_s=100;
%
% Photoreceptor threshold
  u1ort=20;
%
% 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;
%   pause
    options=odeset(options,'JPattern',S)
    [t,u]=ode15s(@pde_1,tout,u0,options);   
  end 
%
% One vector to five 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;           
        if(ncase==1)u1cc(it,ncc)=pcc_s;     end
        if(ncase==2)u1cc(it,ncc)=0.1*pcc_s; end
      end
    end
%
%   u2
    for it=1:nout
      for i=1:nor u2or(it,i)=u(it,nt+i);             end
    end
%
%   Display tabulated numerical solution
%
%   Inner retina
    for it=1:nout
      fprintf('\n    t     zg_ir   u1ir(z,t)')
      for i=1:nir
        fprintf('%5.0f%10.0f%12.2f\n'  ,t(it)/60,zg_ir(i),u1ir(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.0f\n\n',t(it)/60,zg_or(i),u2or(it,i))
      end  
%
%     Fluid layer      
      fprintf('\n    t     zg_fl   u1fl(z,t)\n')
      for i=1:nfl
        fprintf('%5.0f%10.0f%12.2f\n'  ,t(it)/60,zg_fl(i),u1fl(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))
      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');  
%
%   One plot for u2or
    figure(3)
    plot(zg_or,u2or,'-')
    axis([0 200 0 1.0e+07])
    xlabel('zg-or, \mum'); ylabel('u2or(z,t), cells/cm^2');
    title('Outer retina photoreceptors, t=0,20,...,120 min');      
    
    
    

    
    
