% This code simulates invasion percolation ala Kolwankar et al (Europhys
% Lett 62 (4) pp 519-525 (2003))
function hot_vs_cold(H,L,T,tau)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INPUT
% H and L are the height and width of the box
% T is the temperature of the system
% tau is the time you wait before these blocks go away
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
[Y,X]=find(R.^-1 <= tau);
figure(1)
plot(X,Y,'r.');
hold on
[Y2,X2]=find(R.^-1 > tau);
plot(X2,Y2,'b.')
axis equal
axis off