PREPARE_PDE

 create initial data for solving PDE
 
 Input:
     1. ncell
     2. cell_shape
     3. params_domain_geom is a structure with 3 elements:
         Rratio_IN
         include_ECS
         ECS_gap    
     4. params_domain_pde is a structure with 8 elements:
         dcoeff_IN
         dcoeff_OUT
         dcoeff_ECS
         ic_IN
         ic_OUT
         ic_ECS
         kappa_IN_OUT
         kappa_OUT_ECS
 
 Output:
     1. DIFF_cmpts
     2. kappa_bdys
     3. IC_cmpts
     4. OUT_cmpts_index
     5. ECS_cmpts_index
     6. IN_cmpts_index
     7. Ncmpt
     8. Nboundary

0001 function [DIFF_cmpts,kappa_bdys,IC_cmpts,OUT_cmpts_index,ECS_cmpts_index,IN_cmpts_index,Ncmpt,Nboundary] ...
0002     = PREPARE_PDE(ncell,cell_shape,params_domain_geom,params_domain_pde)
0003 
0004 % create initial data for solving PDE
0005 %
0006 % Input:
0007 %     1. ncell
0008 %     2. cell_shape
0009 %     3. params_domain_geom is a structure with 3 elements:
0010 %         Rratio_IN
0011 %         include_ECS
0012 %         ECS_gap
0013 %     4. params_domain_pde is a structure with 8 elements:
0014 %         dcoeff_IN
0015 %         dcoeff_OUT
0016 %         dcoeff_ECS
0017 %         ic_IN
0018 %         ic_OUT
0019 %         ic_ECS
0020 %         kappa_IN_OUT
0021 %         kappa_OUT_ECS
0022 %
0023 % Output:
0024 %     1. DIFF_cmpts
0025 %     2. kappa_bdys
0026 %     3. IC_cmpts
0027 %     4. OUT_cmpts_index
0028 %     5. ECS_cmpts_index
0029 %     6. IN_cmpts_index
0030 %     7. Ncmpt
0031 %     8. Nboundary
0032 
0033 Rratio_IN = params_domain_geom.Rratio_IN;
0034 include_ECS = params_domain_geom.include_ECS;
0035 
0036 if (Rratio_IN <= 0 || Rratio_IN >= 1)
0037     Ncmpt = ncell;
0038 else
0039     Ncmpt = 2*ncell;
0040 end
0041 
0042 if (include_ECS ~= 0)
0043     Ncmpt = Ncmpt + 1;
0044 end
0045 
0046 if (Rratio_IN <= 0 || Rratio_IN >= 1)
0047     if (cell_shape == 2)
0048         % for a cylinder, the side is one boundary, the top and bottom another boundary
0049         Nboundary = 2*ncell;
0050     elseif (cell_shape == 1)
0051         % for a sphere, there is one boundary
0052         Nboundary = ncell;
0053     else
0054         disp('wrong');
0055         stop
0056     end
0057 else
0058     if (cell_shape == 2)
0059         % for a cylinder, there is the axon side wall, the axon top and bottom,
0060         % the myelin side wall, the myelin top and bottom.
0061         Nboundary = 4*ncell;
0062     elseif (cell_shape == 1)
0063         % for a sphere, there is the outer sphere and the inner sphere
0064         Nboundary = 2*ncell;
0065     else
0066         disp('wrong');
0067         stop
0068     end
0069 end
0070 
0071 if (include_ECS ~= 0)
0072     Nboundary = Nboundary+1;
0073 end
0074 
0075 DIFF_cmpts = zeros(1,Ncmpt);
0076 kappa_bdys = zeros(1,Nboundary);
0077 
0078 if (include_ECS ~= 0)
0079     % the ECS is the last cmpt
0080     % the ECS outer boundary is the last boundary
0081     ECS_cmpts_index = Ncmpt;
0082     Box_boundary = Nboundary;
0083 else
0084     ECS_cmpts_index = [];
0085     Box_boundary = [];
0086 end
0087 
0088 DIFF_cmpts(1,ECS_cmpts_index) = params_domain_pde.dcoeff_ECS;
0089 % the outer boundary of ECS is always impermeable
0090 kappa_bdys(1,Box_boundary) = 0;
0091 IC_cmpts(1,ECS_cmpts_index) = params_domain_pde.ic_ECS;
0092 
0093 if (cell_shape == 2)
0094     if (Rratio_IN <= 0 || Rratio_IN >= 1)
0095         OUT_cmpts_index = 1:ncell;
0096         IC_cmpts(1,OUT_cmpts_index) = params_domain_pde.ic_OUT;
0097         DIFF_cmpts(1,OUT_cmpts_index) = params_domain_pde.dcoeff_OUT;
0098         OUT_ECS_boundary = 1:2:2*ncell;
0099         if (include_ECS ~= 0)
0100             kappa_bdys(1,OUT_ECS_boundary) = params_domain_pde.kappa_OUT_ECS;
0101         end 
0102         IN_cmpts_index = [];        
0103     else
0104         IN_cmpts_index = 1:ncell;
0105         IC_cmpts(1,IN_cmpts_index) = params_domain_pde.ic_IN;
0106         DIFF_cmpts(1,IN_cmpts_index) = params_domain_pde.dcoeff_IN;
0107         OUT_cmpts_index = ncell+1:2*ncell;
0108         DIFF_cmpts(1,OUT_cmpts_index) = params_domain_pde.dcoeff_OUT;
0109         IC_cmpts(1,OUT_cmpts_index) = params_domain_pde.ic_OUT;        
0110         IN_OUT_boundary = [1:4:4*ncell];
0111         kappa_bdys(1,IN_OUT_boundary) = params_domain_pde.kappa_IN_OUT;
0112         OUT_ECS_boundary = [3:4:4*ncell];
0113         if (include_ECS ~= 0)
0114             kappa_bdys(1,OUT_ECS_boundary) = params_domain_pde.kappa_OUT_ECS;    
0115         end
0116     end
0117     
0118 elseif (cell_shape == 1)
0119     OUT_cmpts_index = 1:ncell;
0120     IC_cmpts(1,OUT_cmpts_index) = params_domain_pde.ic_OUT;
0121     DIFF_cmpts(1,OUT_cmpts_index) = params_domain_pde.dcoeff_OUT;
0122     if (Rratio_IN <= 0 || Rratio_IN >= 1)        
0123         OUT_ECS_boundary = [1:ncell];
0124         if (include_ECS ~= 0)
0125             kappa_bdys(1,OUT_ECS_boundary) = params_domain_pde.kappa_OUT_ECS;
0126         end
0127         IN_cmpts_index = [];        
0128     else        
0129         IN_cmpts_index = ncell+1:2*ncell;
0130         DIFF_cmpts(1,IN_cmpts_index) = params_domain_pde.dcoeff_IN;
0131         IC_cmpts(1,IN_cmpts_index) = params_domain_pde.ic_IN;        
0132         IN_OUT_boundary = 2:2:2*ncell;
0133         kappa_bdys(1,IN_OUT_boundary) = params_domain_pde.kappa_IN_OUT;
0134         OUT_ECS_boundary = 1:2:2*ncell;
0135         if (include_ECS ~= 0)
0136             kappa_bdys(1,OUT_ECS_boundary) = params_domain_pde.kappa_OUT_ECS;  
0137         end
0138     end
0139 end
0140 
0141 % if (Rratio_IN <= 0)
0142 %     if (cell_shape == 2)
0143 %         % this is the axon
0144 %         IN_ECS_boundary = 1:2:2*ncell;
0145 %         kappa_bdys(1,IN_ECS_boundary) = params_domain_pde.kappa_IN_ECS;
0146 %         OUT_cmpts_index = [];
0147 %     elseif (cell_shape == 1)
0148 %         OUT_ECS_boundary = [1:ncell];
0149 %         kappa_bdys(1,OUT_ECS_boundary) = params_domain_pde.kappa_OUT_ECS;
0150 %         IN_cmpts_index = [];
0151 %     else
0152 %         disp('wrong');
0153 %         stop
0154 %     end
0155 % else
0156 %     if (cell_shape == 2)
0157 %         OUT_cmpt_index = ncell+1:2*ncell;
0158 %         DIFF_cmpts(1,OUT_cmpt_index) = params_domain_pde.dcoeff_OUT;
0159 %         IC_cmpts(1,OUT_cmpt_index) = params_domain_pde.ic_OUT;
0160 %
0161 %         IN_OUT_boundary = [1:4:4*ncell];
0162 %         kappa_bdys(1,IN_OUT_boundary) = params_domain_pde.kappa_IN_OUT;
0163 %         OUT_ECS_boundary = [3:4:4*ncell];
0164 %         kappa_bdys(1,OUT_ECS_boundary) = params_domain_pde.kappa_OUT_ECS;
0165 %
0166 %
0167 %     elseif (cell_shape == 1)
0168 %         IN_cmpts_index = ncell+1:2*ncell;
0169 %         DIFF_cmpts(1,IN_cmpts_index) = params_domain_pde.dcoeff_IN;
0170 %         IC_cmpts(1,IN_cmpts_index) = params_domain_pde.ic_IN;
0171 %
0172 %         IN_OUT_boundary = 2:2:2*ncell;
0173 %         kappa_bdys(1,IN_OUT_boundary) = params_domain_pde.kappa_IN_OUT;
0174 %         OUT_ECS_boundary = 1:2:2*ncell;
0175 %         kappa_bdys(1,OUT_ECS_boundary) = params_domain_pde.kappa_OUT_ECS;
0176 %     else
0177 %         disp('wrong');
0178 %         stop
0179 %     end
0180 % end