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

```