% run and pause tumble process in 1 d 
clear all

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);
cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

N = 1000; % number of particles

% pick some parameter values
kon = 1;% rate of switching to moving
koff = 1; % rate of switching to pause 
v = 1;  % ballistic velocity 
T_max = 100; %maximum time to run simulation
 dt = 10;
 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 = 0; %Initial position
    X = x;
    Xp(1) = x;  % record trajectory for plotting
    Tp(1) = 0;  % record times for plotting
    k = 0; 
    while(t<T_max)
        %for each step we need three 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(2*R(3,1)-1); % direction for the moving state
        t = t + dt1+dt2; %time after swim and tumble
        k=k+1;
        T(k+1) = t;
        x = x+sigma*v*dt1; % swim in direction sigma for time dt1
       
        X(k+1) = x; % save the new position
        if (j==1) % create a trajectory for plotting
        Xp = [Xp,x,x];
        Tp = [Tp,T(k)+dt1,t];
        end
    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;
     
    if (j==1)
        % plot the trajectory
        figure(1)
        plot(Tp,Xp,'linewidth',2)
          xlabel('time','fontsize',20)
        ylabel('x','fontsize',20)
        axis([0 S(end) fix(min(Xp)-1) fix(max(Xp) + 1)])
       if(cmykflg==1)
print('../../figs_c/chapt_4/1_d_run_tmbl_traj','-deps','-cmyk')
 end
    end
    
end
Ys = Ys/N; % is the mean squared displacement (msd)
figure(2)
plot([0,S],[0,Ys],'*',[0,S],[0,S],'--','linewidth',2)
xlabel('time','fontsize',20)
ylabel('Mean Squared Displacement','fontsize',20)       
if(cmykflg==1)
print('../../figs_c/chapt_4/1_d_run_tmbl_msd','-deps','-cmyk')
 end

% 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 = kon*v^2/(koff*(kon+koff));
2*D_est  % D_est is my estimate of the diffusion coefficient
        