% run and tumble process in 2d agent based
clear all
N = 1000; % number of particles

% pick some parameter values
kon = 10;% rate of switching to swimming
koff = 1; % rate of switching to tumbling
rates = [kon,koff];
 
v = 1;% ballistic velocity
T_max = 100; %maximum time to run simulation
 % initial positions
x = zeros(N,1);
y = zeros(N,1);
x1 = x(1);
y1 = y(1);
X = x;  % save the trajectories
Y = y;
 % initial state
 
 % states are:
 % 0 tumbling state
 % 1 moving state
st = ones(N,1); % initially all in the moving state
theta = 2*pi*rand(N,1); %  initial  angle of movement

 dt = 0.1;  % time step 
 t = 0;
 T(1) = 0;
 
Nt = T_max/dt;  % number of simulation time steps 
msd = zeros(Nt,1);
drift = zeros(Nt,2);

for j = 2:Nt  % simulate the motion for each particle
   t = t+dt;
   T(j) = t;
  
    % update the position
    % if in state 1:
    
     x = x+ v.*dt.*st.*cos(theta); % swim in direction theta for time dt
    y = y+ v.*dt.*st.*sin(theta); % swim in direction theta for time dt
    % if in state 0, pick a new angle
     theta = theta+ 2*pi*rand(N,1).*(st==0);
     
    probs = rates(st+1)'*dt;  %These are the probabilities of switching
    st = mod(st+(probs >=rand(N,1)),2);
    
    % if you want, save all the trajectories here, but this step is time
    % consuming
  %   X = [X,x]; % save the new position
   %  Y = [Y,y];
   
   x1 = [x1,x(1)];
   y1 = [y1,y(1)];
     % track the mean square displacement 
     
     msd(j) = sum(x.^2+y.^2)/N;
     drift(j,1) = sum(x)/N;
     drift(j,2) = sum(y)/N;
     
end
   
        % plot a sample trajectory
        figure(1)
        plot( x1,y1,'linewidth',2 )
    title('Sample trajectory','fontsize',20)
    xlabel('x','fontsize',20)
    ylabel('y','fontsize',20)
    
    figure(2)
    plot(T, msd,'*')
    xlabel('time','fontsize',20)
    ylabel('MSD','fontsize',20)
    
    
    figure(3)
    % final positions
    plot(x,y,'*')
    title('Final Positions','fontsize',20)
    xlabel('x','fontsize',20)
    ylabel('y','fontsize',20)
    figure(4)
    plot(drift(:,1),drift(:,2))
    title('Mean Displacement','fontsize',20)
    xlabel('x','fontsize',20)
    ylabel('y','fontsize',20)
       