function tlr2_tumor_code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
% This is a companion code for the results shown in Fig.5B "Analysis  %
% of ligand molecule movement through the tissue space  and its role  %
% in receptor saturation" from a paper: Karolak et al. "Targeting     % 
% Ligand Specificity Linked to Tumor Tissue Topological Heterogeneity % 
% via Single-Cell Micro-Pharmacological Modeling", Scientific Reports %
% 2018, 8:3638                                                        %
%                                                                     %
% Author: Aleksandra Karolak, Moffitt Cancer Center                   %
%                                                                     %
% output files: a directory "Output" is created and the following     %
%   text files will be saved:                                         %
%      mol_in_#.txt -- with coordinates of all internalized molecules %
%      mol_4ins_#.txt -- with the number of a cell (1 column) and the %
%                 coordinates of molecules internalized by this cell  %
%      mol_4out_#.txt -- with coordinates of all particles that left  %
%                 the domain                                          %
%                                                                     %
% adjustable parameters:                                              %
%   1. Particles release scheme (3 hrs for slow and 1 min for fast)   %
%   2. Concentration of released biomarker (fenestration at time)     %
%   3. Diffusion coefficient (pp in the code)                         %
%   4. Binding affinity (frequency of binding, probability)           %
%                                                                     %
% An example from Figure 5B requires the following parameters:        %
%   releaseScheme=1;  fenestr=10;  pp=4;  bindingAffinity=10;         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% RELEASE
  releaseScheme   = 1;            % 1 for fast, 0 for slow
    if releaseScheme == 1           % FAST release
      fenestrationIter = 1; 
      finalFenIter     = 500;
    elseif releaseScheme == 0       % SLOW release
      fenestrationIter = 250; 
      finalFenIter     = 125000;
    end

% CONCENTRATION
  fenestr = 10;           % 1 for conc 500, 2 for 1000 and so on (x500)

% DIFFUSION
  pp = 4;                 % 4 for diff 10-6;   3 for 10-5;   2 for 10-4.

% AFFINITY 
  bindingAffinity = 10;   % 10 for aff 100%, 100 for 10%, 1000 for aff 1%.

% TIME       
  NmaxIter    = 216000;   % xxx/600 for minutes
  Nsave       = 500;      % frequency of saving 
  Ndraw       = 250;      % frequency of drawing images
  cell_count  = 8;        % number of cells in the tissue
  pathOut=['Output']; mkdir(pathOut)  % output directory for saving data
       
% DOMAIN PARAMETERS
  xmin=0; xmax=0.21; ymin=0; ymax=0.21;                   % boundaries
  [ygr,xgr] = meshgrid(xmin:0.002:xmax,ymin:0.002:ymax);  
  Nx=size(xgr,1); Ny=size(xgr,2);
  xx=xgr(:,1);    yy=ygr(1,:);    hg=(xgr(2,1)-xgr(1,1));

  maxCellLoad        = 4;        % if zero, only one ligand will bind;
  distToReceptor     = 0.0055;   % for the binding to take place  
  distToCellBoundary = 0.00575;  % for particles from inside to outside
  rad                = 0.0099;  
  dt                 = 0.1;     
  fenY=[ymin:((ymax-ymin)/(fenestr-1)):ymax]; fenX=xmin*ones(size(fenY));    
 
%  PARAMETERS CONTROLLING THE CELLULAR UPTAKE TIME 
   time4uptake = 1500;
   randLimit   = 0.999;
   flowLimit   = 0.75;          % When to start uptake
 
   maxPart=fenestr; maxInfl=10*10; 
   diff=[2.5e-3,2.5e-4,2.5e-5,2.5e-6,2.5e-7,2.5e-8];
   permDist=0.120; 
   peds=10000;
   dif=diff(pp); 

  % --- READ COORDINATES OF EACH CELL:
   for k=1:cell_count
     xy_cells{k}=DefineCells(k);
   end

  % --- CALCULATE NUMBER OF RECEPTORS AND APPROXIMATE THE CELL CENTER 
   for cell_number = 1:cell_count                  %cell_count
     nn=0;
     Ncells(cell_number) = size(xy_cells{cell_number},1);
     cellLoad{cell_number} = zeros(Ncells(cell_number),1);
     centers(cell_number,:) = [sum(xy_cells{cell_number}...
             (nn+1:nn+Ncells(cell_number),1)),...
	 sum(xy_cells{cell_number}(nn+1:nn+Ncells(cell_number),2))]/Ncells(cell_number);
     nn=nn+Ncells(cell_number);
   end

   xpart=zeros(maxInfl*maxPart,1); ypart=zeros(maxInfl*maxPart,1); 
   numPart=0;  
   xins=zeros(1,4); yins=zeros(1,4);
   x4ins=zeros(1,2); y4ins=zeros(1,2);
   x4out=zeros(1,1); x4out=zeros(1,2);
   NNi=0; Nthr=0; Fthr=0; thrPart=0; iter = -1; thrPart=0; 

   if (Ndraw>0)
     fig1=figure('position',[200,200,500,500])
     fig2=figure('position',[800,200,500,500])
   end
      
% --- MAIN LOOP OVER ITERATIONS:
   while ((iter<NmaxIter))  
     if mod(iter,50)==0
       [iter,NmaxIter]  
     end
       
     iter=iter+1;                                      
     if ((mod(iter,fenestrationIter)==0) && iter < finalFenIter)
       xpart(numPart+1:numPart+maxPart,1)=fenX;
       if (fenestr == 1)
         ypart(numPart+1:numPart+maxPart,1)=(fenY+0.21)*rand;  % only for fenestration = 1
       elseif (fenestr == 2)
         fenY(1)=fenY(2);
         ypart(numPart+1:numPart+maxPart,1)=fenY*rand;  % only for fenestration = 2
       elseif (fenestr == 5)
	     fenY(1:4)=fenY(5);
         ypart(numPart+1:numPart+maxPart,1)=fenY*rand;  % only for fenestration = 5
       else
         ypart(numPart+1:numPart+maxPart,1)=fenY;
       end
       numPart=numPart+maxPart;
     end
     xtemp=xpart(1:numPart,1); ytemp=ypart(1:numPart,1);
     xmemin=xins; ymemin=yins;

     % --- CALCULATE FLOW RATE:
     indbefore =find(xpart < 0.1);
     indafter  =find(xpart >=0.1);
     flow=length(indbefore)/length(indafter);

     for cell_number=1:cell_count
       for iji=1:Ncells(cell_number)    % ***
         clear xtemp ytemp  
         xtemp=xpart(1:numPart,1);
         ytemp=ypart(1:numPart,1);

         dx = sqrt((xpart(1:numPart,1)-xy_cells{cell_number}(iji,1)).^2+...
                   (ypart(1:numPart,1)-xy_cells{cell_number}(iji,2)).^2);
         indfind =find(dx<distToReceptor);

         if length(indfind)>0
           dxmin = min(dx); 
           if ((mod(iter,bindingAffinity)==0) && ...
               (cellLoad{cell_number}(iji,1) <= maxCellLoad))
             ind1 = find(dx == dxmin);
             ind2 = find(dx ~= dxmin);
             inn=length(ind1);
             outn=length(ind2);
             xins(NNi+1:NNi+inn,1)=xpart(ind1,1);  
             yins(NNi+1:NNi+inn,1)=ypart(ind1,1);  
             for ll=1:inn  
	       xins(NNi+ll,2)=xy_cells{cell_number}(iji,1);
               yins(NNi+ll,2)=xy_cells{cell_number}(iji,2);
               x4ins(NNi+ll,2)=xy_cells{cell_number}(iji,1);
	       y4ins(NNi+ll,2)=xy_cells{cell_number}(iji,2);
               x4ins(NNi+ll,1)=cell_number;
	       y4ins(NNi+ll,1)=cell_number;
             end
             cellLoad{cell_number}(iji,1)=cellLoad{cell_number}(iji,1)+1;
             NNi=NNi+inn;

             clear xpart ypart;
             xpart=xtemp(ind2); 
             ypart=ytemp(ind2);
             numPart=size(xpart,1);
           end 
         end
       end    % iji=1:Ncells(cell_number) ***
         
       for j = 1:NNi
         dcent = sqrt((xins(j,2)-centers(cell_number,1)).^2 +...
                      (yins(j,2)-centers(cell_number,2)).^2);
         indin = find(dcent<=3*rad);
         if ((length(indin)>0) && (xins(j,4) == 0) && ...
            (xins(j,3) > time4uptake) && (rand > randLimit) )
		   xins(j,2)=centers(cell_number,1)+...
                     0.5*rand*((xins(j,2)-centers(cell_number,1))); 
           yins(j,2)=centers(cell_number,2)+...
                     0.5*rand*((yins(j,2)-centers(cell_number,2))); 
           xins(j,4)=1;                       
         end  
       end  
     end  % cell_number=1:cell_count

     xins(:,3)=xins(:,3)+1*dt;                % aging of particle 

     % --- diffusion component 
     wkx=sqrt(2*dif*dt)*randn(1,peds);      
     wky=sqrt(2*dif*dt)*randn(1,peds); 
     selx=round((peds-1)*rand(1,numel(xpart(1:numPart,1)))+1); 
     sely=round((peds-1)*rand(1,numel(ypart(1:numPart,1)))+1); 
     difx=wkx(selx)';        
     dify=wky(sely)';
  
     % --- particle movement: memorize old position and add diffusion component
     nn=0;
     xmem=xpart; ymem=ypart; 
     xpart(1:numPart,1)=xpart(1:numPart,1)+difx; 
     ypart(1:numPart,1)=ypart(1:numPart,1)+dify;
     iinn=0;

     % --- Crossing cell boundaries
     % --- verify if particles move too close to cell or cross cell boundary
     for cell_number=1:cell_count
       for i=1:Ncells(cell_number)
         distb=sqrt((xpart(1:numPart,1)-xy_cells{cell_number}(i,1)).^2+...
                    (ypart(1:numPart,1)-xy_cells{cell_number}(i,2)).^2);
         indcp=find(distb<distToCellBoundary);

         if length(indcp)>0     %--- if molecule is found within distance from cell boundary
   		   for k=1:length(indcp)
             if i == 1
               distbck(k)=sqrt((xpart(indcp(k),1)-xy_cells{cell_number}...
                  (Ncells(cell_number),1)).^2+(ypart(indcp(k),1)-...
                  xy_cells{cell_number}(Ncells(cell_number),2)).^2);
               distfrw(k)=sqrt((xpart(indcp(k),1)-xy_cells{cell_number}...
                  (i+1,1)).^2+(ypart(indcp(k),1)-xy_cells{cell_number}(i+1,2)).^2);

               if distbck(k) < distfrw(k)
 			     v1x(k) = xy_cells{cell_number}(i,1)-xy_cells{cell_number}...
                     (Ncells(cell_number),1);
                 v1y(k) = xy_cells{cell_number}(i,2)-xy_cells{cell_number}...
                     (Ncells(cell_number),2);
                 v2x(k) = xpart(indcp(k),1)-xy_cells{cell_number}(i,1);
                 v2y(k) = ypart(indcp(k),1)-xy_cells{cell_number}(i,2);
                 crossP(k) = (v1x(k) * v2y(k) - v1y(k) * v2x(k));
		       else
 			     v1x(k) = xy_cells{cell_number}(i+1,1)-xy_cells{cell_number}(i,1);
                 v1y(k) = xy_cells{cell_number}(i+1,2)-xy_cells{cell_number}(i,2);
                 v2x(k) = xpart(indcp(k),1)-xy_cells{cell_number}(i+1,1);
                 v2y(k) = ypart(indcp(k),1)-xy_cells{cell_number}(i+1,2);
                 crossP(k) = -(v1x(k) * v2y(k) - v1y(k) * v2x(k));
               end  
                
             elseif (i > 1) && (i < Ncells(cell_number))	           
               distbck(k)=sqrt((xpart(indcp(k),1)-xy_cells{cell_number}(i-1,1)).^2+...
			       (ypart(indcp(k),1)-xy_cells{cell_number}(i-1,2)).^2);
               distfrw(k)=sqrt((xpart(indcp(k),1)-xy_cells{cell_number}(i+1,1)).^2+...
                   (ypart(indcp(k),1)-xy_cells{cell_number}(i+1,2)).^2);

               if distbck(k) < distfrw(k)
 			     v1x(k) = xy_cells{cell_number}(i,1)-xy_cells{cell_number}(i-1,1);
                 v1y(k) = xy_cells{cell_number}(i,2)-xy_cells{cell_number}(i-1,2);
                 v2x(k) = xpart(indcp(k),1)-xy_cells{cell_number}(i,1);
                 v2y(k) = ypart(indcp(k),1)-xy_cells{cell_number}(i,2);                       
                 crossP(k) = (v1x(k) * v2y(k) - v1y(k) * v2x(k));
               else
  			     v1x(k) = xy_cells{cell_number}(i+1,1)-xy_cells{cell_number}(i,1);
                 v1y(k) = xy_cells{cell_number}(i+1,2)-xy_cells{cell_number}(i,2);
                 v2x(k) = xpart(indcp(k),1)-xy_cells{cell_number}(i+1,1);
                 v2y(k) = ypart(indcp(k),1)-xy_cells{cell_number}(i+1,2);
                 crossP(k) = -(v1x(k) * v2y(k) - v1y(k) * v2x(k));  
               end
                
             else     %if i == Ncells(cell_number)
               distbck(k)=sqrt((xpart(indcp(k),1)-xy_cells{cell_number}(i-1,1)).^2+...
				  (ypart(indcp(k),1)-xy_cells{cell_number}(i-1,2)).^2);
               distfrw(k)=sqrt((xpart(indcp(k),1)-xy_cells{cell_number}(1,1)).^2+...
                  (ypart(indcp(k),1)-xy_cells{cell_number}(1,2)).^2);

               if distbck(k) < distfrw(k)
 			     v1x(k) = xy_cells{cell_number}(i,1)-xy_cells{cell_number}(i-1,1);
                 v1y(k) = xy_cells{cell_number}(i,2)-xy_cells{cell_number}(i-1,2);
                 v2x(k) = xpart(indcp(k),1)-xy_cells{cell_number}(i,1);
                 v2y(k) = ypart(indcp(k),1)-xy_cells{cell_number}(i,2);  
                 crossP(k) = (v1x(k) * v2y(k) - v1y(k) * v2x(k));
               else
                 v1x(k) = xy_cells{cell_number}(1,1)-xy_cells{cell_number}(i,1);
                 v1y(k) = xy_cells{cell_number}(1,2)-xy_cells{cell_number}(i,2);
                 v2x(k) = xpart(indcp(k),1)-xy_cells{cell_number}(1,1);
                 v2y(k) = ypart(indcp(k),1)-xy_cells{cell_number}(1,2);
                 crossP(k) = -(v1x(k) * v2y(k) - v1y(k) * v2x(k));                    
               end  
             end
              
             if (crossP(k) <= 0) % < for inside the cell 
               x4out(indcp(k),1)=xmem(indcp(k),1); x4out(indcp(k),2)=ymem(indcp(k),1);
               th=2*pi*rand; rr=5*rad;
               xpart(indcp(k),1)= xmem(indcp(k),1);
               ypart(indcp(k),1)= ymem(indcp(k),1);  
             end
           end  % for k=
         end   % length     
       end  % i 

	   dist=sqrt((xpart(1:numPart,1)-centers(cell_number,1)).^2+...
                 (ypart(1:numPart,1)-centers(cell_number,2)).^2);
       indx=find(dist<1.5*rad);
                                    
       if size(indx,1)>0       
         xpart(indx,1)=xmem(indx,1);       
         ypart(indx,1)=ymem(indx,1);    
       end
          
     end   % cell count

     clear xmem ymem        
     ind1=(ypart(1:numPart,1)<ymin); if (sum(ind1)>0); ypart(ind1)=ymin+0.001; end             
     ind1=(ypart(1:numPart,1)>ymax); if (sum(ind1)>0); ypart(ind1)=ymax-0.001; end             
     ind1=(xpart(1:numPart,1)<0);    if (sum(ind1)>0); xpart(ind1)=xmin+0.001; end 
     thrPart=sum(xpart>=permDist);
     if (thrPart>=0.25*maxPart*maxInfl)&&(Fthr==0)
       Nthr=iter; Fthr=1;  
     end
 
      
     
     
     % -- draw the figure with rendered drug concentration 
     if (mod(iter,Ndraw)==0)
       xy=[xpart(1:numPart,1),ypart(1:numPart,1)];
       xyins=[xins(1:NNi,1),yins(1:NNi,1),xins(1:NNi,2),yins(1:NNi,2)];
       Nin=size(xyins,1);

       
       figure(fig1);
         clf
         colordef black
             
         conc=zeros(Nx,Ny);
         for kk=1:size(xy,1)
           if ((xy(kk,1)>=xmin)&&(xy(kk,1)<=xmax)&&...
               (xy(kk,2)>=ymin)&&(xy(kk,2)<=ymax))  
             idx=1+floor((xy(kk,1)-xmin)/hg);
             idy=1+floor((xy(kk,2)-ymin)/hg);
             conc(idx,idy)=conc(idx,idy)+1;
           end
         end
      
         for kk=1:Nin
           if ((xyins(kk,3)>=xmin)&&(xyins(kk,3)<=xmax)&&...
	           (xyins(kk,4)>=ymin)&&(xyins(kk,4)<=ymax))  
	         idx=1+floor((xyins(kk,3)-xmin)/hg);
             idy=1+floor((xyins(kk,4)-ymin)/hg);
             conc(idx,idy)=conc(idx,idy)+1;  
           end
         end       
 
         pcolor(xx,yy,conc')
         %imagesc(xx,yy,conc')
         hold on;
         view(2)
         axis equal
         axis([xmin,xmax,ymin,ymax])
         colorbar
         colormap(jet)
         title(['TLR2 concentration in the in vivo tissue; ',...
                '  Iter=',num2str(iter)],'FontSize',15)
         shading interp
         caxis([0,5])   
      
         set(gcf,'renderer','painters')
         for kk=1:8
           plot(xy_cells{kk}(:,1),xy_cells{kk}(:,2),'.','Color',...
               [1,1,0],'MarkerSize', 0.35)
         end
         pause(0.1)




        % figure without rendering
          figure(fig2);
          clf
          colordef white

	  fill([xmin,xmax,xmax,xmin,xmin],[ymin,ymin,ymax,ymax,ymin],'w')
	  axis([xmin,xmax,ymin,ymax])
          axis equal
	  axis([xmin,xmax,ymin,ymax])
          hold on

	  for kk=1:8
		   plot(xy_cells{kk}(:,1),xy_cells{kk}(:,2),'.','Color',...
			[0,0,0],'MarkerSize', 0.35)
          end

          plot(xy(:,1),xy(:,2),'r.','MarkerSize',15)
          title(['Tissue in vivo ','  Iter=',num2str(iter)],'FontSize',15)
          if (Nin>0)
              plot(xyins(:,3),xyins(:,4),'g.','MarkerSize',15)
          end
          pause(0.01)

        
        % SAVING THE FIGURES
        if Nsave>0
          figure(fig1)  
          filename=[pathOut,'/Fig_',num2str(iter)];
          print('-djpeg100','-painters','-r300',filename)       
          
          figure(fig2)  
          filename=[pathOut,'/Fig2_',num2str(iter)];
          print('-djpeg100','-painters','-r300',filename)       
        end
      end
      
      
      if (mod(iter,Nsave)==0)   
	    datSave=[xpart(1:numPart,1),ypart(1:numPart,1)];                 
        filenam=['mol_',num2str(iter),'.txt'];
        pathname=[pathOut,'/',filenam];
        save(pathname,'datSave','-ascii','-double')
        clear datSave
          
        datSave=[xins(1:NNi,1),yins(1:NNi,1),xins(1:NNi,2),yins(1:NNi,2)];
        filenam=['mol_in_',num2str(iter),'.txt'];
        pathname=[pathOut,'/',filenam];
        save(pathname,'datSave','-ascii','-double')
        clear datSave

    	forout = find(x4out(:,1) > 0);
        datSave=[x4out(forout,1),x4out(forout,2)];
        filenam=['mol_4out_',num2str(iter),'.txt'];
        pathname=[pathOut,'/',filenam];
        save(pathname,'datSave','-ascii','-double')
        clear datSave
        
        forin = find(x4ins(:,1) > 0);
        datSave=[x4ins(forin,1),x4ins(forin,2),y4ins(forin,2)];
        filenam=['mol_4ins_',num2str(iter),'.txt'];
        pathname=[pathOut,'/',filenam];
        save(pathname,'datSave','-ascii','-double')
        clear datSave  
      end  
    end   % while
      

    clear centers xy_cells bpts fpts xpart ypart xins yins

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DefineCells -- define coordinates of all cells
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cell=DefineCells(CellNum)
  switch CellNum
    case 1,
    cell=[0.110791, 0.0407772; 0.112677, 0.0414060; 0.113935, 0.0420349;...
          0.116241, 0.0432926; 0.117918, 0.0447600; 0.119385, 0.0468562;...
          0.120014, 0.0487428; 0.120433, 0.0506294; 0.120643, 0.0523063;...
          0.120853, 0.0539833; 0.120853, 0.0554507; 0.121691, 0.0579661;...
          0.122110, 0.0611104; 0.122320, 0.0640451; 0.122529, 0.066351;...
          0.123368, 0.0682376; 0.123787, 0.0715915; 0.122739, 0.0751551;...
          0.120853, 0.0776705; 0.119804, 0.0795571; 0.118966, 0.0814437;...
          0.117289, 0.0835399; 0.114564, 0.0856361; 0.111, 0.0860554;...
          0.107437, 0.0860554; 0.105131, 0.0860554; 0.101148, 0.0860554;...
          0.099262, 0.0854265; 0.0967461, 0.0850073; 0.0952787, 0.0839592;...
          0.0925536, 0.0820726;0.0910863, 0.080186; 0.0900382, 0.078509;...
          0.0896189, 0.0764128;0.0894093, 0.074107; 0.0887805, 0.0713819;...
          0.087942, 0.0701242;0.0885708, 0.0671895; 0.0896189, 0.0644644;...
          0.0896189, 0.0615297; 0.0904574, 0.0583854; 0.0912959, 0.0550314;...
          0.092344, 0.0514679; 0.0936017, 0.0483235; 0.0961172, 0.0449696;...
          0.097794, 0.0432926; 0.099681, 0.0418253; 0.101987, 0.0409868;...
          0.104921, 0.0403579; 0.108066, 0.0403579];
     case 2,
     cell=[0.142458, 0.0662319; 0.144295, 0.0656809; 0.145948, 0.0654972;...
          0.148887, 0.0653135; 0.151458, 0.0653135; 0.154213, 0.0656809;...
          0.155683, 0.0660483; 0.156601, 0.066783; 0.158989, 0.0682524;...
          0.160642, 0.0697218; 0.16303, 0.0713749; 0.164316, 0.073579;...
          0.165969, 0.0755995; 0.165969, 0.0779873; 0.16505, 0.0800077;...
          0.165234, 0.0822118; 0.164132, 0.0845996; 0.162846, 0.0864364;... 
          0.161377, 0.0886405; 0.15954, 0.0897426; 0.156601, 0.0913957;...
          0.154581, 0.0919467; 0.152009, 0.0935998; 0.149805, 0.0950692;... 
          0.144662, 0.0981917; 0.142642, 0.0987428; 0.140437, 0.0991101;... 
          0.133458, 0.0991101; 0.12905, 0.0985591; 0.124825, 0.0987428;... 
          0.121519, 0.097457; 0.120784, 0.0943345; 0.120784, 0.0904773;... 
          0.122253, 0.0880895; 0.124458, 0.0860691; 0.127396, 0.0827629;... 
          0.129233, 0.0796404; 0.129233, 0.0776199; 0.131621, 0.0754158;... 
          0.133825, 0.073579; 0.135295, 0.0719259; 0.136397, 0.0704565;... 
          0.137499, 0.0693544; 0.138784, 0.067885; 0.141356, 0.0665993];     
     case 3,
     cell=[0.178678, 0.0877944; 0.184412, 0.0875896; 0.190351, 0.0888184;...
          0.193218, 0.0920952; 0.195266, 0.0978296; 0.194652, 0.10295;...
          0.193014, 0.106841; 0.190146, 0.111756; 0.187894, 0.114214;...
          0.187894, 0.116671; 0.18605, 0.11831; 0.182774, 0.120972;...
          0.180111, 0.123225; 0.178678, 0.125478; 0.17622, 0.12773;...
          0.173967, 0.128959; 0.170486, 0.130188; 0.167209, 0.130188;...
          0.162089, 0.129369; 0.156354, 0.12855; 0.153897, 0.127526;...
          0.151439, 0.125682; 0.149186, 0.123839; 0.148572, 0.12261;...
          0.148162, 0.120767; 0.147958, 0.118924; 0.147958, 0.116876;...
          0.147958, 0.113599; 0.148777, 0.110937; 0.149801, 0.10725;...
          0.153078, 0.104178; 0.155126, 0.10213; 0.156969, 0.100082;...
          0.158198, 0.0980344; 0.160246, 0.0972152; 0.161679, 0.0957816;...
          0.163113, 0.0947576; 0.164546, 0.0937336; 0.16598, 0.0925048;...
          0.168028, 0.0916856; 0.169257, 0.091276; 0.1711, 0.090252;...
          0.172738, 0.0900472; 0.173762, 0.0896376; 0.175606, 0.0888184;...
          0.177244, 0.088204];     
     case 4,
     cell=[0.148062, 0.13781; 0.150513, 0.139816; 0.15341, 0.141822;...
          0.155639, 0.143382; 0.156307, 0.146279; 0.156976, 0.149844;...
          0.155639, 0.154524; 0.153187, 0.156976; 0.14873, 0.159873;...
          0.144942, 0.16121; 0.140262, 0.162993; 0.13625, 0.162993;...
          0.130234, 0.162547;0.127114, 0.161878; 0.124439, 0.15965;...
          0.120428, 0.156976; 0.117754, 0.154301; 0.116417, 0.152741;...
          0.114188, 0.150959; 0.11196, 0.149844; 0.109063, 0.147616;...
          0.107948, 0.145833; 0.107948, 0.141376;0.107726, 0.13781;...
          0.107726, 0.136028; 0.108171, 0.133353; 0.109063, 0.131348;...
          0.109063, 0.129342; 0.110623, 0.127336; 0.111291, 0.125108;...
          0.112851, 0.123325;0.115303, 0.122657; 0.119314, 0.122434;...
          0.123102, 0.122657; 0.125331, 0.123102; 0.128005, 0.124439;...
          0.129788, 0.125554; 0.132016, 0.126891; 0.134022, 0.128896;...
          0.136028, 0.130234; 0.138479, 0.132462; 0.140708, 0.133576;...
          0.143159, 0.134913; 0.14717, 0.136473];     
     case 5,
     cell=[0.0531326, 0.172141; 0.0543274, 0.173336; 0.0564782, 0.175009;...
          0.05839, 0.177638; 0.0591069, 0.180505; 0.0591069, 0.18409;...
          0.0588679, 0.186719; 0.0562392, 0.190064; 0.0533715, 0.191737;...
          0.049548, 0.192693; 0.0462023, 0.192693; 0.0428567, 0.192693;...
          0.0397501, 0.192693; 0.0356875, 0.192454; 0.0318639, 0.191976;...
          0.0289963, 0.189825; 0.0268455, 0.187674; 0.0251727, 0.183851;...
          0.0251727, 0.180505; 0.0249337, 0.178593; 0.0249337, 0.176921;...
          0.0256506, 0.174053; 0.0270845, 0.172619; 0.0280404, 0.170229;...
          0.0306691, 0.168318; 0.0335368, 0.167362; 0.0356875, 0.166406;...
          0.0385552, 0.166406; 0.0411839, 0.166884; 0.0430957, 0.166884;...
          0.0442905, 0.167362; 0.0459634, 0.168079; 0.0466803, 0.168557;...
          0.0483531, 0.169512; 0.0497869, 0.16999; 0.0512208, 0.171185];     
     case 6,
     cell=[0.0361784, 0.0704722; 0.0369004, 0.0728787; 0.0369004, 0.0762479;...
          0.0369004, 0.0810611; 0.036419, 0.084671; 0.0373817, 0.0890028;...
          0.0388256, 0.0921314; 0.0405102, 0.094538; 0.0409915, 0.0974259;...
          0.0407509, 0.100073; 0.0400289, 0.10248; 0.0381036, 0.103683;...
          0.0356971, 0.105849; 0.0323278, 0.106571; 0.0304026, 0.106571;...
          0.0277553, 0.106571; 0.0253488, 0.10609; 0.0229422, 0.105127;...
          0.0205356, 0.102239; 0.018851, 0.0995918; 0.0171664, 0.0969445;...
          0.0152411, 0.094538; 0.0150005, 0.0918907; 0.0145192, 0.0887622;...
          0.0145192, 0.0865962; 0.0147598, 0.083949; 0.0150005, 0.0815424;...
          0.0154818, 0.0788952; 0.0162038, 0.0772106; 0.0169257, 0.074804;...
          0.0186103, 0.0726381; 0.0200543, 0.0707128; 0.0212576, 0.0692689;...
          0.0222202, 0.0678249; 0.0236642, 0.0671029; 0.0248674, 0.066381;...
          0.0260707, 0.0658996; 0.0275147, 0.065659; 0.028718, 0.0654183;...
          0.0296806, 0.0654183; 0.0313652, 0.0654183; 0.0323278, 0.065659;...
          0.0337718, 0.0658996; 0.0344938, 0.0668623; 0.0352157, 0.0673436;...
          0.0356971, 0.0680656];     
     case 7,
     cell=[0.00667467, 0.0367042; 0.0058379, 0.0350306; 0.0058379, 0.0327295;...
          0.0058379, 0.03001; 0.0058379, 0.0266629; 0.00730225, 0.0249894;...
          0.00730225, 0.0218515; 0.00918498, 0.019132; 0.0108585, 0.0176676;...
          0.0137872, 0.0176676; 0.0167159, 0.0174584; 0.0194354, 0.0176676;...
          0.0219458, 0.019132; 0.0244561, 0.0208055; 0.0271756, 0.0228974;...
          0.0288491, 0.0262445; 0.0294767, 0.028964; 0.0294767, 0.0314744;...
          0.0305227, 0.0337755; 0.0305227, 0.0362858; 0.0301043, 0.0398421;...
          0.0296859, 0.0410972; 0.0282216, 0.0440259; 0.0278032, 0.0459087;...
          0.0267572, 0.0469546; 0.0240377, 0.0475822; 0.0219458, 0.0475822;...
          0.0196446, 0.0477914; 0.0183895, 0.0475822; 0.0173435, 0.047373;...
          0.0158792, 0.0471638; 0.0144148, 0.0471638; 0.013578, 0.0467455;...
          0.0121137, 0.0463271; 0.0106493, 0.0456995; 0.0100218, 0.0444443;...
          0.00960337, 0.0431892; 0.0087666, 0.0417248; 0.00813902, 0.0404697;...
          0.00772063, 0.0385869];     
     case 8,
     cell=[0.0619857, 0.0221274; 0.064362, 0.0229916; 0.0673865, 0.0253679;...
          0.0697629, 0.0279603; 0.0719232, 0.0307688; 0.0730034, 0.0344414;...
          0.0742996, 0.03833; 0.0727873, 0.0415705; 0.0714911, 0.0441629;...
          0.0660903, 0.0454591; 0.0613376, 0.0458911; 0.0563688, 0.0461072;...
          0.0520481, 0.0461072; 0.0490236, 0.0458911; 0.0464312, 0.0443789;...
          0.0412464, 0.0389781; 0.0401662, 0.0355215; 0.0395181, 0.032497;...
          0.0395181, 0.0290405; 0.0397342, 0.0266641; 0.0401662, 0.0247198;...
          0.0405983, 0.0225595; 0.0421105, 0.0210473; 0.0436228, 0.019103;...
          0.0459992, 0.0184549; 0.0479435, 0.0184549; 0.0503198, 0.0184549;...
          0.0522641, 0.0188869; 0.0561527, 0.019319; 0.0574489, 0.0206152;...
          0.0583131, 0.0206152];         
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%