% run and tumble process in 1 d 

% this code is set up to simulate many different cases of run and tumble in
% 1 d
% choose parameters to fit the case you desire:
clear all
N = 10000; % number of particles

% pick some parameter values
kon = 1;% rate of switching to moving
koff = 1; % rate of switching to pause 
Vm = 1;% ballistic velocity to the left
Vp = 1;  % ballistic velocity to the right
V = [Vm,Vp];
alpha = 1/2 ; % 0< alpha <1 probability of going left, 1-alpha = probability of going right
% alpha = 1/2 -> equal probability of going left or right
D = 0; %diffusion coefficient for motion while paused
T_max = 100; %maximum time to run simulation
 dt = 10; %time step for sampling
 ns = fix(T_max/dt); % number of time steps
Ys =zeros(1,ns); %store the mean squared displacements
Ym = zeros(1,ns); % store the mean displacement
for j = 1:N  % simulate the motion for each particle
    clear X
    
    T(1) = 0;  %initial time
    t = 0;
    x = 0; %Initial position
    X = x;
     
    k = 0; 
    while(t<T_max)
        %for each step we need four random numbers
        R = rand(3,1); 
        dt1 = -log(R(1,1))/koff; % amount of time in moving state
        dt2 = -log(R(2,1))/kon; % amount of time in paused state
        sigma = sign(R(3,1)-alpha); % direction for the moving state
        sdx = (sigma +3)/2;
        Rn = randn(1,1); 
        t = t + dt1+dt2; %time after swim and diffuse
        k=k+1;
        T(k+1) = t;
        x = x+sigma*V(sdx)*dt1 +sqrt(2*D*dt2)*Rn; % swim in direction sigma for time dt1
       % and diffuse for time dt2
        X(k+1) = 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));
      
     %Estimate the position by linear interpolation (although not exactly
     %correct)
     % save the msd
     Ys = Ys+Xs.^2;
     % save the mean
     Ym = Ym + Xs;
     
    %in order to determine the msd
    
    % not done yet.
    
    
    if (j==1)
        % plot the trajectory
        figure(1)
        plot(T,X,'linewidth',2)
          xlabel('time','fontsize',20)
        ylabel('x','fontsize',20)
        
        
        
    end
    
end
Ym = Ym/N;  % this is the mean;
Ys = Ys/N; % is the mean squared displacement (msd)
Ys = Ys-Ym.^2; %this is the variance
    
     

% find the slope of the msd
p = polyfit(S,Ys,1);
p(1) %is the slope of the msd = 2 D_eff
% compare this with an estimate for D
%D_est = v^2*kon*(4*alpha*(1-alpha)*kon*(kon+2*koff)+koff^2)/(kon+koff)^3/koff + D*koff/(kon+koff);
D_est = ((Vp+Vm)^2*kon^2*alpha*(1-alpha)*(kon + 2*koff)+koff^2*kon*(Vm^2*alpha+Vp^2*(1-alpha))) ...
    /(koff*(kon+koff)^3) + D*koff/(kon+koff);
2*D_est  % D_est is my estimate of the diffusion coefficient

Mv = -(Vm*alpha-Vp*(1-alpha))*kon/(koff+kon);  % estimate of mean velocity
figure(2)
plot([0,S],[0,Ym],'*',[0,S],Mv*[0,S],'--')
xlabel('time','fontsize',20)
ylabel('Mean Displacement','fontsize',20)  

figure(3)
plot([0,S],[0,Ys],'*',[0,S],2*D_est*[0,S],'--')
xlabel('time','fontsize',20)
ylabel('Variance','fontsize',20)   