Fortran Program Code for 2D Elastic and Inelastic Collision of 2 Balls

Copy and Paste from Browser

      
c******************************************************************************
c   This program is a 'remote' 2D-collision detector for two balls on linear
c   trajectories and returns, if applicable, the location of the collision for 
c   both balls as well as the new velocity vectors (assuming a partially elastic
c   collision as defined by the restitution coefficient).
c
c   All variables apart from 'error' are of Double Precision Floating Point type.
c   In 'free' mode no positions but only the initial velocities and an impact 
c   angle are required.
c   All variables apart from 'mode' and 'error' are of Double Precision
c   Floating Point type.
c
c   The Parameters are:
c
c    mode  (character*4) (if='free' alpha must be supplied; otherwise arbitrary)
c    alpha (impact angle) only required in mode='free'; 
c                       should be between -PI/2 and PI/2 (0 = head-on collision))
c    R    (restitution coefficient)  between 0 and 1 (1=perfectly elastic collision)
c    m1    (mass of ball 1)
c    m2    (mass of ball 2)
c    r1    (radius of ball 1)       not needed for 'free' mode
c    r2    (radius of ball 2)                "
c  * x1    (x-coordinate of ball 1)          "
c  * y1    (y-coordinate of ball 1)          "
c  * x2    (x-coordinate of ball 2)          "
c  * y2    (y-coordinate of ball 2)          "
c  * vx1   (velocity x-component of ball 1) 
c  * vy1   (velocity y-component of ball 1)         
c  * vx2   (velocity x-component of ball 2)         
c  * vy2   (velocity y-component of ball 2)
c  * error (Integer*2) (0: no error 
c                       1: balls do not collide 
c                       2: initial positions impossible (balls overlap)
c
c   In the calling program, the arguments corresponding to the parameters
c   with an asterisk (*) will be updated (the positions and velocities however
c   only if 'error'=0).
c   All variables should have the same data types in the calling program
c   and all should be initialized before calling the subroutine.
c
c   This program is free to use for everybody. However, you use it at your own
c   risk and I do not accept any liability resulting from incorrect behaviour.
c   I have tested the program for numerous cases and I could not see anything 
c   wrong with it but I can not guarantee that it is bug-free under any 
c   circumstances.
c   I do in general not recommend the use of single precision because for
c   iterative computations the rounding errors can accumulate and render the
c   result useless if there are not enough decimal places; however for uncritical
c   applications or if speed is of paramount importance, you can change the 
c   program to single precision by changing the variable declaration to REAL*4
c   and remove the 'D' in front of the mathematical function names.
c
c
c   I would appreciate if you could report any problems to me
c   (for contact details see  http://www.plasmaphysics.org.uk/feedback.htm ).
c
c   Thomas Smid, January  2004
c                December 2005 (corrected faulty collision detection; 
c                               a few minor changes to improve speed;
c                               added simplified code without collision detection)
c                December 2009 (generalization to partially inelastic collisions)
c**********************************************************************************
                       
     
       SUBROUTINE collision2D(mode,alpha,R,m1,m2,r1,r2,
     &   x1,y1,x2,y2,vx1,vy1,vx2,vy2,error)
   
       IMPLICIT REAL*8 (A-Z)
       CHARACTER*4 mode
       INTEGER*2 error
       

c      **** initialize some additional variables ****
       pi2=2.0D0*DACOS(-1.0D0);
       error=0
       r12=r1+r2
       m21=m2/m1
       x21=x2-x1
       y21=y2-y1
       vx21=vx2-vx1
       vy21=vy2-vy1

       vx_cm = (m1*vx1+m2*vx2)/(m1+m2) 
       vy_cm = (m1*vy1+m2*vy2)/(m1+m2) 

c      **** return if relative velocity =0  ****
       IF (vx21 .EQ. 0 .AND. vy21 .EQ. 0) THEN
       error=1
       RETURN
       END IF
       
c      *** calculate relative velocity angle             
       gammav=DATAN2(-vy21,-vx21)



c******** this block only if initial positions are given ***

       IF (mode .NE. 'free') THEN


       
       d=DSQRT(x21*x21 +y21*y21)
       
c      **** return if distance between balls smaller than sum of radii ***
       IF (d .LT. r12) THEN
       error=2
       RETURN
       END IF
       
c      **** calculate relative position angle and normalized impact parameter ***
       gammaxy=DATAN2(y21,x21)
       dgamma=gammaxy-gammav 
         IF (dgamma .GT. pi2) THEN
         dgamma=dgamma-pi2
           ELSE 
           IF (dgamma .LT. -pi2) THEN
           dgamma=dgamma+pi2
           END IF
         END IF  
       dr=d*DSIN(dgamma)/r12
       
       
c      **** return old positions and velocities if balls do not collide ***
       IF (  (dgamma .GT. pi2/4.0D0  .AND.  ddgamma .LT. 0.75D0*pi2  ) 
     &         .OR.  DABS(dr) .GT. 1.0D0 ) THEN
       error=1
       RETURN
       END IF
c      **** calculate impact angle if balls do collide ***              
       alpha=DASIN(dr)
       

       
c      **** calculate time to collision ***      
       dc=d*DCOS(dgamma)
       IF (dc .GT. 0) THEN
       sqs=1.0D0
       ELSE
       sqs=-1.0D0
       END IF      
       t=(dc-sqs*r12*DSQRT(1-dr*dr))/ DSQRT(vx21*vx21+vy21*vy21)
     
c      **** update positions *** 
       x1=x1+vx1*t
       y1=y1+vy1*t
       x2=x2+vx2*t
       y2=y2+vy2*t
       


       END IF
       
c******** END 'this block only if initial positions are given' ***
      
       
       
c  ***  update velocities ***

       a=DTAN(gammav+alpha)

       dvx2=-2*(vx21+a*vy21)/((1+a*a)*(1+m21))
       
       vx2=vx2+dvx2
       vy2=vy2+a*dvx2
       vx1=vx1-m21*dvx2
       vy1=vy1-a*m21*dvx2       
       

c  ***  velocity correction for inelastic collisions ***
	   
       vx1=(vx1-vx_cm)*R + vx_cm
       vy1=(vy1-vy_cm)*R + vy_cm
       vx2=(vx2-vx_cm)*R + vx_cm
       vy2=(vy2-vy_cm)*R + vy_cm
       
       
       RETURN
       END





c******************************************************************************
c  Simplified Version
c  The advantage of the 'remote' collision detection in the program above is 
c that one does not have to continuously track the balls to detect a collision. 
c  The program needs only to be called once for any two balls unless their 
c velocity changes. However, if somebody wants to use a separate collision 
c  detection routine for whatever reason, below is a simplified version of the 
c  code which just calculates the new velocities, assuming the balls are already 
c  touching (this condition is important as otherwise the results will be incorrect)
c****************************************************************************
     
       SUBROUTINE collision2Ds(R,m1,m2,x1,y1,x2,y2,vx1,vy1,vx2,vy2)
   
       IMPLICIT REAL*8 (A-Z)
       CHARACTER*4 mode
       INTEGER*2 error
       
     
       m21=m2/m1
       x21=x2-x1
       y21=y2-y1;
       vx21=vx2-vx1
       vy21=vy2-vy1

       vx_cm = (m1*vx1+m2*vx2)/(m1+m2) 
       vy_cm = (m1*vy1+m2*vy2)/(m1+m2) 



c      *** return old velocities if balls are not approaching ***
       IF ( (vx21*x21 + vy21*y21) .GE. 0) THEN
               RETURN
       END IF


c      *** I have inserted the following statements to avoid a zero divide; 
c      *** (for single precision calculations, 
c      ***         1.0D-12 should be replaced by a larger value). **********  
  
       fy21=1.0D-12*DABS(y21)                            
       IF ( DABS(x21) .LT. fy21 ) THEN
                  IF (x21 .LT. 0) THEN 
                       sign=-1.0D0  
                  ELSE  
                       sign=1.0D0
                  END IF      
                  x21=fy21*sign
        END IF 


       
       
c  ***  update velocities ***

       a=y21/x21

       dvx2=-2*(vx21 +a*vy21)/((1+a*a)*(1+m21)) 
       
       vx2=vx2+dvx2
       vy2=vy2+a*dvx2
       vx1=vx1-m21*dvx2
       vy1=vy1-a*m21*dvx2       
       
c  ***  velocity correction for inelastic collisions ***
	   
       vx1=(vx1-vx_cm)*R + vx_cm
       vy1=(vy1-vy_cm)*R + vy_cm
       vx2=(vx2-vx_cm)*R + vx_cm
       vy2=(vy2-vy_cm)*R + vy_cm        
       
       RETURN
       END




C++ Code

(Back To) Description