% direction switching 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
kswitchm  = 1;% rate of switching  from right to left
kswitchp  = 1;% rate of switching  from left to right
% for this code, these are assumed to be equal
v = 1;  % ballistic velocity 
T_max = 100; %maximum time to run simulation
 dt = 10;  % time increment for sampling the process
 ns = fix(T_max/dt); % number of time steps
ms =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 two random numbers
        R = rand(2,1); 
        dt1 = -log(R(1,1))/kswitchm; % amount of time   moving right
        dt2 = -log(R(2,1))/kswitchp; % amount of time  moving left
        
        t = t + dt1+dt2; %time after two switches
        k=k+1;
        T(k+1) = t;
        x = x+v*dt1-v*dt2; % move to the right for time dt1 and to the left for time dt2
       
        X(k+1) = x; % save the new position
        if (j==1) % create a trajectory for plotting
        Xp = [Xp,X(k)+v*dt1,x];
        Tp = [Tp,T(k)+dt1,T(k)+dt1+dt2];
        end
    end
    % now sample this trajectory at equal time steps
   
    S = [1:ns]*dt; % times at which trajectory is sampled
    for jj = 1:ns
        Sj(jj) = min(find(S(jj)<T)); % find the points that surround the desired time value
    end
     
     J = [1:ns];
     %Estimate the position by linear interpolation  
     Xs=X(Sj-1)+(S(J)-T(Sj-1)).*(X(Sj)-X(Sj-1))./(T(Sj)-T(Sj-1));
      % save this trajectory
      if (j==1)
          Xsv = Xs;
      else
          Xsv = [Xsv;Xs];
      end
  
   
    if (j==1)
        % plot the trajectory
        figure(1)
        plot(Tp,Xp, 'linewidth',2)
        xlabel('time','fontsize',20)
        ylabel('x','fontsize',20)
        title('Sample Trajectory','fontsize',18)
        axis([0 S(end) fix(min(Xp)-1) fix(max(Xp) + 1)])
       % title('Sample Trajectory','fontsize',20)
         if(cmykflg==1)
print('../../figs_c/chapt_4/dir_switch_traj','-deps','-cmyk')
 end
    end
    
end

% Now post processing on Xsv to find msd
msd = sum(Xsv.^2)/N % is the mean squared displacement  
 
% find the slope of the msd  (fit to a line through the origin)
R = sum(msd.*S);
X = sum(S.^2);
sl = R/X;

sl %is the slope of the msd = 2 D_eff
% compare this with an estimate for D
D_est = v^2/(2*kswitchm);
2*D_est  % D_est is my estimate of the diffusion coefficient

figure(2)
plot([0,S],[0,msd],'*', [0,S],sl*[0,S],'--','linewidth', 2)
 xlabel('time','fontsize',20)
        ylabel('Mean Squared Displacement','fontsize',20)   
  if(cmykflg==1)
print('../../figs_c/chapt_4/dir_switch_msd','-deps','-cmyk')
 end      