% run and tumble process in 2d agent based with chemoattractant
function run_tumble_chemo
 global Lx Ly
 
 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

cmykflg=1;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

N = 1000; % number of particles
 Lx = 5;  % location of the target
 Ly = 3;
% pick some parameter values
kon = 10;% rate of switching to swimming
koff = 1; % rate of switching to tumbling
rates = [kon,koff];
chilist = [0,0.1];   % do two simulations; one with oure diffusuon, one with chemoattraction.
for j = 1:2
chi = chilist(j); %strength of chemoattractant   
 v = 1;% ballistic velocity
T_max = 50; %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.05;  % 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:
    c_old = external(x,y);
     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
    c_new = external(x,y);
    c_grad = (c_new-c_old)/dt;
    % if in state 0, pick a new angle
     theta = theta+ 2*pi*rand(N,1).*(st==0);
     
    probs = (rates(st+1)'-chi*c_grad)*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 
     
    
     drift(j,1) = sum(x)/N;
     drift(j,2) = sum(y)/N;
      msd(j) = sum((x-drift(j,1)).^2+(y-drift(j,2)).^2)/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)
    hold on
    figure(2)
    plot(T, msd,'*')
    xlabel('time','fontsize',20)
    ylabel('Variance','fontsize',20)
    hold on
    if(cmykflg==1)
  print('../../figs_c/chapt_11/chemotax_agent_variance','-deps','-cmyk')
     end 
    
    figure(3)
    % final positions
    plot(x,y,'*')
    title('Final Positions','fontsize',20)
    xlabel('x','fontsize',20)
    ylabel('y','fontsize',20)
    hold on
    if(cmykflg==1)
  print('../../figs_c/chapt_11/chemotax_agent_positions','-deps','-cmyk')
    end 
 
    figure(4)
    plot(drift(:,1),drift(:,2),'linewidth',2)
    title('Mean Displacement','fontsize',20)
    xlabel('x','fontsize',20)
    ylabel('y','fontsize',20)
     hold on
     
 
end

   
    function c = external(x,y)
    % 
    global Lx Ly
    c =  -((x-Lx).^2+(y-Ly).^2) ;
    
       