% run and tumble process in 2d uses Gillespie algorithm
clear all
N = 5000; % number of particles

% pick some parameter values
kon = 1;% rate of switching to swimming
koff = 1; % rate of switching to tumbling
D = 1;  %Diffusion coefficient for tumble
v = 1;% ballistic velocity
T_max = 1000; %maximum time to run simulation
 dt = 10; % time step to sample data
 ns = fix(T_max/dt); % number of time steps
Ys =zeros(1,ns); %store the mean squared displacements
for j = 1:N  % simulate the motion for each particle
    clear X
    
   T(1) = 0;  %initial time
    t = 0;
    x = zeros(2,1); %Initial position in 2 d
    theta = 2*pi*rand(1,1); %initial angle
    X = x;
   
    k = 0; 
    while(t<T_max)
        %for each step we need three random numbers
        R = rand(2,1); % uniformly distributed
        %Rn = randn(1,1); % normally distributed
        dt1 = -log(R(1,1))/koff; % amount of time in the moving state
        dt2 = -log(R(2,1))/kon; % amount of time in the tumbling state
        
        t = t + dt1+dt2; %time after swim and tumble
        k=k+1;
        T(k+1) = t;
        x = x+ v*dt1*[ cos(theta);sin(theta)]; % swim in direction theta for time dt1
       %theta = theta + sqrt(2*D*dt2)*Rn; % if theta diffuses for time dt2
       theta = 2*pi*rand(1,1);  %if the new angle is uniformly distributed
        X = [X,x]; % save the new position
         
        end
  
    % now sample this trajectory at equal time steps
   
    S = [1:ns]*dt;
    for jj = 1:ns
        Sj(jj) = min(find(S(jj)<T));
    end
     
     J = [1:ns];
     %Xs=X(Sj-1)+(S(J)-T(Sj-1)).*(X(Sj)-X(Sj-1))./(T(Sj)-T(Sj-1));
     Xs = X(:,Sj-1) ;
     %Estimate the position by linear interpolation (although not exactly
     %correct)
     % save the msd
     Ys = Ys+sum(Xs.^2);
    
    if (j==1)
        % plot the trajectory
        figure(1)
        plot(X(1,:),X(2,:),Xs(1,:),Xs(2,:),'*')
        
    end
end
    
    figure(2)
    plot(S, Ys/N,'*')
    
    p=polyfit(S,Ys/N,1);
    p(1)  % this is the slope of the msd
        