% This code simulates invasion percolation ala Kolwankar et al (Europhys
% Lett 62 (4) pp 519-525 (2003))

function [t,N_bin]=invasion_perc(H,L,T,bins)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% INPUT
% H and L are the height and width of the box

% T is the temperature of the system 

% We are interested in the question: "how many sites dissolved between time
% t and time t+delta_t?". the histogram of this data can be normalized to
% tell us the probability that a N-avalanche will occure (i.e. N sites will
% dissolve in a time delta_t. "bins" is defined such that 
% delta_t=(total run time)/bins

% OUTPUT

% t is a vector such that t(i) is the time at which the ith site dissolved

% N_bins is a vector such that N_bins(i) is the number of sites that
% dissolved between delta_t*i and delta_t*(i+1). the histogram of this
% vector can be used to construct the pdf for the avalanche sizes

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if T<.001
    error('The temperature is to low for this code to work. Increase temperature to 10^(-3).')
end




% initialize the system
R=ones(H,L+1);

R(:,2:(end-1))=exp(-rand(H,L-1)/T); % give random reaction rates to all of the sites
R(:,1)=zeros(H,1); % the sites that have been removed
R(:,end)=zeros(H,1);
front_x=2*ones(H,1);% the x_coordinates of the front
front_y=(1:H)';% the x_coordinates of the front

dead_x=ones(H,1);% the x_coordinates where further material cannot be removed
dead_y=[(1:H)';(1:H)']; %the x_coordinates where further material cannot be removed
dead_x=[dead_x;(L+1)*ones(H,1)];
% run the simulation until half of the particles are gone
T=floor(H*(L-1)/2)+1;
dead_x=[dead_x;ones(T,1)];
dead_y=[dead_y;ones(T,1)];


%% run the simulation until you get to the far wall
%kill=0; % kill will change from zero to one when a point on the far wall gets hit
ind=0;
t=zeros(T,1);
figure
for k=2:T
    dead_k=(k-1)+2*H;
%% Find the amount of time until the next cell goes
I=sub2ind([H,L],front_y,front_x);
r=R(I);
delta_t=sum(r)^-1;
t(k)=t(k-1)+delta_t;
    
%% find the site the disappeared

% Select a point at random such that the probability of selecting it is
% R_i/sum(R) where R_i is its reactivity
    p=rand(1)/delta_t;% pick a random numebr between 0 and sum(r)
    c=cumsum(r);
    % find the point in reaction rates that it is nearest to
    pt=c-p;
    index=sum(pt<0); % this is the index of the largest element smaller than p
    index=index+1;
    [ROW,COL]=ind2sub([H,L+1],I(index));

    % find the points around the one that will be removed
    
    neigh_row=[ROW+1;ROW-1;ROW;ROW];
    neigh_col=[COL;COL;COL+1;COL-1];
    kill_this=find(neigh_row<1 | neigh_row>H);
    neigh_col(kill_this)=[];% take out the points that are out of bounds
    neigh_row(kill_this)=[];
    pos_i=sub2ind([H,L+1],neigh_row,neigh_col);
   
    % remove the point
    front_x(index)=[];
    front_y(index)=[];
    dead_x(dead_k)=COL;
    dead_y(dead_k)=ROW;
    dead_i=sub2ind([H,L+1],dead_y(1:dead_k),dead_x(1:dead_k));
    useless_i=[dead_i;I];% these indicies are either dead or on the front

    for j=1:length(pos_i)
              
        % add this point if it is not already dead or on the front        
        if sum(pos_i(j)==useless_i)==0;
            front_y=[front_y;neigh_row(j)];
            front_x=[front_x;neigh_col(j)];
        end
    end
    
%     [Y,X]=find(-log(R)<10 & R>0);
%     [Y_b,X_b]=find(-log(R)>10 & R>0);
    
    if floor(ind/10)==ind/10
       plot(dead_x,dead_y,'.k')
%         hold on
%         plot(X_b,Y_b,'.b')
%         hold off
%         axis equal
      axis([0,L,0,H])
    drawnow
    end
    ind=ind+1;
end
% at this point I have taken out ind many sites

%% divide the total time into 1000 segments
t_bin=linspace(t(end)/10,t(end),bins+1);% cut off the early stages in the evolution
start=t_bin(1:(end-1));
fin=t_bin(2:end);

N_bin=ones(1,bins);
for i=1:bins
    b_1=start(i);
    b_f=fin(i);
    N_bin(i)=sum(b_1<t & t<=b_f);
end