% This code simulates a simple poisson jump process
clear all

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

N = 1000; % number of particles

T_max = 1000; %maximum time to run the 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
Xmn = zeros(1,ns); % store the mean value
for j = 1:N  % simulate the motion for each particle, one at a time
    clear X T
    
   T(1) = 0;  %initial time
    t = 0;
    x = 0; %Initial position  
    
    X = x;
   
    k = 0; 
    while(t<T_max)
        %for each step we need two uniformly distributed random numbers
        R = rand(2,1); % uniformly distributed
       
        delt = - log(R(1,1)); %this is the time step
        r = R(2,1); 
        dx = (-log(1-r)*(r>1/2)+log(r)*(r <=1/2)); % the distance moved in this step
        % or
        %dx = 1;  %to make it deterministic
        t = t + delt; %time after this step
        k=k+1;
        T(k+1) = t; 
        x = x+ dx; % move x by the amount dx 
        dxs(k+1) = dx;
        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) ;
      Xmn = Xmn+Xs;  %accumulate the mean 
%      % save the msd
      Ys = Ys+ Xs.^2;
%     
      if (j<=10)
        % plot  sample trajectories
        figure(1)
        plot(T,X )
        hold on
      end
    
end
    Vs  = Ys/N-Xmn.^2/N^2; % the variance
 figure(2)
plot(S,Vs,'*')   
    p=polyfit(S,Vs,1);
    p(1)  % this is the slope of the msd
% 
figure(3)
plot(S,Xmn/N,'*')