%%% Simple model for lateral inhibition %%%



%%%%%%%%%%% Input parameters %%%%%%%%%%%%%%

g = 0.1;          %%%%% Nonlinear dissipation
eps = -2.0;       %%%%% Coupling strength
Omega = 0.0;     %%%%% Production rate
delta = 0.5;      %%%%% Topology parameter (delta = 0.0 square lattice) (delta = 0.5 hexagonal lattice)
N = 10000;        %%%%% Number of time steps


%%%%%%%%%%%%%%%%%%%%%%%%% Do not change from here %%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%% Integration parameters %%%%%%%%%%%

dt = 0.01;          %%%%%%% Time step
j = 9;              %%%%%%% This parameter determines the size of the lattice
size = 3 + 2*(j-1); %%%%%%% Actual size of lattice
center = j + 1;     %%%%%%% Center of lattice
counter = 0;        %%%%%%% Initializing counter for movie


%%%%%%%% Setting up initiual conditions %%%%%%%%%%  

 [X,Y] = meshgrid(1:size,1:size);    %%%%%%% Setting mesh
 u  =  ones(size,size);              %%%%%%% Initial condition
 Input = Omega*ones(size,size);      %%%%%%% Setting the input
    

    
%%%%%%%%%%%%%%%% Numerical Integration %%%%%%%%%%%%%%%%%%%%
    
    for i = 1 : N
        
        
 
 
        %%%%% Corner I %%%%%%
        
        u(1,1) = (-1*u(1,1)-g*u(1,1)^3+eps*(u(1,2)+u(2,1))+Input(1,1))*dt + u(1,1);
        
        %%%%% Corner II %%%%%%
        
        u(1,end) = (-1*u(1,end)-g*u(1,end)^3+eps*(u(1,end-1)+u(2,end))+Input(1,end))*dt + u(1,end);
        
        %%%%% Corner III %%%%%%
        
        u(end,1) = (-1*u(end,1)-g*u(end,1)^3+eps*u(end-1,1)+Input(end,1))*dt + u(end,1);
        
        u(end,1) = u(end,1) + 2*delta*eps*dt*u(end-1,2);
        
        %%%%% Corner IV %%%%%%
        
        u(end,end) = (-1*u(end,end)-g*u(end,end)^3+eps*u(end-1,end-1)+Input(end,end))*dt + u(end,end);
        
        u(end,end) =  u(end,end) + 2*delta*eps*dt*u(end-1,end-1);
        
        %%%%% Line I %%%%%%
        
        u(1,2:end-1) = (-1*u(1,2:end-1)-g*u(1,2:end-1).^3+eps*(u(1,1:end-2)+u(1,3:end)+u(2,2:end-1))+Input(1,2:end-1))*dt + u(1,2:end-1);
        
        u(1,2:end-1) = u(1,2:end-1)+ delta/2*eps*dt*(u(2,1:end-2)+u(2,3:end));
        
        %%%%% Line II %%%%%%
        
        u(2:end-1,1) = (-1*u(2:end-1,1)-g*u(2:end-1,1).^3+eps*(u(1:end-2,1)+u(3:end,1)+u(2:end-1,2))+Input(2:end-1,1))*dt + u(2:end-1,1);
        
        u(2:end-1,1) =  u(2:end-1,1) + 2*delta*eps*dt*u(1:end-2,2);
           
        %%%%% Line III %%%%%
        
        u(end,3:2:end-2) = (-1*u(end,3:2:end-2)-g*u(end,3:2:end-2).^3+eps*(u(end-1,3:2:end-2))+Input(end,3:2:end-2))*dt + u(end,3:2:end-2);
      
        u(end,3:2:end-2) =  u(end,3:2:end-2) + 2*delta*eps*dt*(u(end-1,2:2:end-3)+u(end-1,4:2:end-1));
        
        %%%%% Line IV %%%%%%
        
        u(2:end-1,end) = (-1*u(2:end-1,end)-g*u(2:end-1,end).^3+eps*(u(1:end-2,end)+u(3:end,end))+Input(2:end-1,end))*dt + u(2:end-1,end);
         
        u(2:end-1,end) =  u(2:end-1,end) + 2*delta*eps*dt*u(1:end-2,end-1);
        
        
        %%%%%%% Bulk %%%%%%%
        
         
        u(2:end-1,2:end-1) = (-1*u(2:end-1,2:end-1)-g*u(2:end-1,2:end-1).^3 + eps*(u(1:end-2,2:end-1)+u(3:end,2:end-1)+u(2:end-1,1:end-2)+u(2:end-1,3:end))+Input(2:end-1,2:end-1))*dt + u(2:end-1,2:end-1);
        
        u(2:end-1,2:2:end-1) = u(2:end-1,2:2:end-1) + 2*delta*eps*dt*(u(3:end,1:2:end-2) + u(3:end,3:2:end));
        
        u(2:end-1,3:2:end-2) = u(2:end-1,3:2:end-2)+ 2*delta*eps*dt*(u(1:end-2,2:2:end-3) + u(1:end-2,4:2:end-1));
        
        %%%%%% Time evolution %%%%%
        
        U_d(i,1) = u(2,center); %center
        
        
        if i ==1 
            
            counter = counter +1;
       
       Frame_U(counter,:,:) = u; 
       
        else if round(i/10)==i/10
       
       counter = counter +1;
       
       Frame_U(counter,:,:) = u;
       
            end
       
          
        end
     
   
      
          
   end
     
      
        

  %%%%%%%%%%%%%%%% Setting the topology of the meshgrid %%%%%%%%%%%%%%%  
    
  [X,Y] = meshgrid(1:size,1:size); 
  Y(:,2:2:end) = delta + Y(:,2:2:end);                       %%%%% Shift to hexagonal lattice
  x = X(:);                            
  y = Y(:);
  U = u(:);
    
  [v,c] = voronoin(horzcat(x,y));                            %%%%% Creating Voronoi diagram
  mymap_red = horzcat((0:.01:1)',zeros(101,1),zeros(101,1)); %%%%% Setting up color scale
 

  

 %%%%%%%%%%%% Plotting and creating movie %%%%%%%%%%%%% 
 
 
 
figure
for j = max(counter)
    % Example of plot
    
    W = squeeze(Frame_U(j,:,:));
   
   for i = 1:length(c) 
if all(c{i}~=1)   % If at least one of the indices is 1, 
                  % then it is an open region and we can't 
                  % patch that.
                             
patch(v(c{i},1),v(c{i},2),W(i)); % use color i.
end
end
  colormap(mymap_red);
 caxis([min(U) max(U)])
 colorbar
 xlim([2.0 size-1.0])
 ylim([2.0 size-1.0])
 set(gca,'Yticklabel',[]) 
 set(gca,'Xticklabel',[])
  set(gca,'Color',[1, 1, 1])
    

   
end


