function [ v1, v2 ] = collide(  u1, u2 , im1, im2  ) 
  ## returns the new velocities of two particles 
  ## with inverse-masses im1 and im2 
  ## after an elastic collision. 
  ##
  ## Dimensions: u1,u2:   [L/T]
  ##             im1,im2: [1/M]
  ##             v1,v2:   [L/T]

  ## compute centre of mass velocity (see proof below)
  c         = ( im2 * u1 + im1 * u2  ) / ( im1 + im2 ) ;
  ##  [L/T] = ( [L/T/M]  +  [L/T/M]  ) / ( [1/M]+[1/M] )  

  ## in the CM frame, velocities reverse
  v1       = c - ( u1 - c ) ;
  v2       = c - ( u2 - c ) ;
  ## [L/T] = [L/T] - [L/T] + [L/T] 
  
endfunction

## See CHECKcollide for a script that checks various special cases.

## centre of mass velocity:
## m_TOT c = m1 u1 + m2 u2
##  c = ( m1 u1 + m2 u2 ) / ( m1 + m2 )
##    =  (im2 u1 + im1 u2) / (m1+m2)/(m1m2)
##    =  (im2 u1 + im1 u2) / (im1 + im2 ) 

