% 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