## Given a state (x,v) of length N and inverse-masses im(1..N)
## where x(1) = left-hand wall  and im(1) = 0 
## and   x(N) = right hand wall and im(N) = 0 

## determine which of the inter-particle collisions will take place next
## -- or is the next event one of the data-collection times {Tc, 2 Tc, 3Tc...}? --
## and simulate the motion until that event, and implement the collision. 

time_to_next_clock = Tc ;

j = 0 ; ## counter: will run up to Number_samples_required 
logbook = zeros(Number_samples_required , 2*N ) ; 

########################################################################
## set up graph
figure(1);
y = zeros(size(x)) ;
## look at the current velocities to guess a good range for the velocity graph
vmax=max(v);
vmin=max(-v) ;
vmax=max(vmax,vmin);
axis([x(1)-1 x(N)+1 -vmax*5.8 vmax*5.8]); 
gset border 15 
gset pointsize 3
########################################################################

################################################################################
## Main loop: 
##             Keep simulating the motion until we have all the data we want
while ( j < Number_samples_required ) 

  ## run through all (left-hand) particles finding the time until collision with their 
  ## right neighbours.
  ## Initialize the collisiontime vector to large positive numbers. 
  collisiontime=ones(1,N) * 2 * time_to_next_clock ;
  ## The Nth particle does not have a right neighbour, so we can use its "collision time" 
  ## to store the time of the next clock event, which is when we will gather data for histograms
  collisiontime(1,N) =  time_to_next_clock ; ## this is used to indicate a clock event
  for n=1:N-1
    dv = v(1,n) - v(1,n+1) ;
    ## if the right neighbour is coming towards us (in our frame of reference)
    ## then work out when the collision would be
    if ( dv > 0.0 ) 
      collisiontime(n) = (x(1,n+1)-x(1,n))/dv ; 
    endif
  endfor

  ## decide which collision needs to be rendered next
  [a,b] = min( collisiontime ) ;

  ## simulate the motion for a time given by the min
  dt = a ; 
  x = x + v * dt ;
  time_to_next_clock = time_to_next_clock - dt ;
  
  left = b ;  ## b holds the index of the particle that has its collision next
  if( left < N ) 
    right = b+1;
    ## change the velocities of the dudes involved in a collision
    [ v(1,left), v(1,right) ] = collide(  v(1,left), v(1,right) , im(1,left), im(1,right) ) ;
##    plot( x , y , "x6" , x , v , "o3");
  else
    ## no collision, it's clock time, so time to take a sample!
    j ++ ;
    if(doplot) 
      plot( x , y , "x5" , x , v , "o3");
      pause(0.075);
    endif
    time_to_next_clock = Tc ;
    logbook(j,:) = [x,v ] ; 
  endif
endwhile

