clear all;
  
  %Generation of bacteria run data
  theta=zeros(300,1);                   %Declares 300x1 matrice for turn angle
  theta=randraw('vonmises',[0,1],300);  %Generates tumbling angle
  
  x=zeros(300,1);
  y=zeros(300,1);
  
  for n=1:300
      plot(x,y,'b.-');
      v_run(n)=normrnd(10,1);            %Generates run speed
      t_run(n)=normrnd(1,0.1);           %Generates run duration
      t_tumb(n)=normrnd(0.1,0.01);       %Generates tumbling duration    
      M(n)=getframe;
      if n==1;
          beta(n)=0;                     %Initialise cummulative angle to 0
          x(n)=0;                        %Initialise origin to (x,y)=(0,0)
          y(n)=0;                   
      else
          beta(n)=beta(n-1)+theta(n-1);  %Generates cummulative turn angle
          x(n)=x(n-1)+v_run(n)*cos(beta(n))*t_run(n);    %Displacement=Velocity*Time
          y(n)=y(n-1)+v_run(n)*sin(beta(n))*t_run(n);   
      end;  
  end;
  
  movie2avi(M,'Motility.avi');
  
  figure
  subplot(2,3,1); hist(v_run,300);  %Plots histogram of run velocity
  Title('PDF of Run Velocities');
  xlabel('Velocity / (um/s)'); ylabel('Frequency');
  
  subplot(2,3,2); hist(t_run,300);  %Plots histogram of run duration
  Title('PDF of Run Duration');
  xlabel('Time / s'); ylabel('Frequency');
  
  subplot(2,3,3); hist(theta,300);  %Plots histogram of turn angle
  Title('PDF of Tumbling Angle');
  xlabel('Tumbling Angle'); ylabel('Frequency');
    
  subplot(2,3,4); hist(t_tumb,300); %Plots histogram of tumbling duration
  Title('PDF of Run Duration');
  xlabel('Time / s'); ylabel('Frequency');
  
  subplot(2,3,5); plot(x,y);        %Plots 2D displacement
  Title('Movement Trace by Single Bacteria');
  xlabel('x-displacement / um');
  ylabel('y-displacement / um');
  
  %Calculate statistical data from synthetic data generated
  
  v_run_mu=mean(v_run)      %Calculate mean run velocity
  t_run_mu=mean(t_run)      %Calculate mean run duration
  t_tumb_mu=mean(t_tumb)    %Calculate mean tumbling time
  theta_mu=mean(theta)      %Calculate mean turn angle
  
  v_run_sigma=std(v_run)    %Calculate run velocity standard deviation
  t_run_sigma=std(t_run)    %Calculate run duration standard deviation
  t_tumb_sigma=std(t_tumb)  %Calculate tumbling duration standard deviation
  theta_sigma=std(theta)    %Calculate turn angle standard deviation