% an agent based simulation of SIRS epidemic with diffusion 
clear
% set some parameters
L = 10; % length and width of 2D domain
dt = 0.01; %the time step  
t = 0.
Nstep = 200;  % the number of time steps to take
D = 1;  % the diffusion coefficient for the particles
scD = sqrt(2*D*dt);  % variance for diffusion
%there are three possible reactions
rxn = 3;  % reaction (transmission) rate I+S -> 2I
% set rxn = 0 to see pure diffusion
Rxn = 1; %size of the reaction neighborhood
rmv = 0.5; % recovery rate I -> R
%rmv = 0; % for no recovery
% set rvm = 0 to see Fisher's equation in an agent based model
rs = 0.; % R -> S rate 
% set rs >0 to see S->I->R->S disease dynamics

 %assign colors to the different particle types
col(1) = 'g';
col(2) = 'r';
col(3) = 'b';

p = 1.0; %This is the density of the initial susceptible population

% this means that there are a total of p*L^2 particles 
N = fix(p*L^2); % total number of particles
dist = zeros(N,N); % the interparticle distances
% randomly assign initial positions
x = rand(N,2)*L;
%there are three possible states
% s=1; susceptible
% s=2 ; infected
% s=3; recovered
s = ones(N,1); %initially all are susceptible

% add one infected at the center of the domain
x(N,:) = [L/2,L/2];
s(N) = 2;

for n = 1:Nstep

% plot the initial configuration
figure(1)
for k = 1:N
    plot(x(k,1),x(k,2), strcat(col(s(k)),'*'),'linewidth',2);
    hold on
end
 
hold off

%there are two steps to take
%1) move  all the particles by diffusion

x = x + scD*randn(N,2);
x = abs(x); % reflecting boundaries at x=0;
x = L - abs(x-L); % reflecting boundaries at x=L;

% calculate reaction probabilities

for j = 1:N
% reaction probability for susceptibles:    
    if(s(j)==1)
        % calculate distances between susceptible and infected particles
        dist = zeros(1,N);
    for k = 1:N
        if (s(k) ==2)
        dist(k) = sum((x(j,:)-x(k,:)).^2);
        end
    end
    DltR = sum((dist < Rxn^2)&(dist>0));
    rxp(j) = rxn*DltR/Rxn^2;
    end
    if (s(j)==2)
        rxp(j) = rmv;
    end
    if (s(j) ==3)
        rxp(j) = rs;
    end
end
 
rxpp = find(rxp*dt>rand(1,N)); %These are the particles that change state
s(rxpp) = s(rxpp) + 1;
rxpp = find(s==4);
s(rxpp) = 1;
% one cycle of the process is complete
end






