Fortran Program Code for 2D Elastic and Inelastic Collision of 2 Balls
Copy and Paste from Browser
c******************************************************************************
c This program is a 'remote' 2Dcollision 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 = headon 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 (xcoordinate of ball 1) "
c * y1 (ycoordinate of ball 1) "
c * x2 (xcoordinate of ball 2) "
c * y2 (ycoordinate of ball 2) "
c * vx1 (velocity xcomponent of ball 1)
c * vy1 (velocity ycomponent of ball 1)
c * vx2 (velocity xcomponent of ball 2)
c * vy2 (velocity ycomponent 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 bugfree 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 (AZ)
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=x2x1
y21=y2y1
vx21=vx2vx1
vy21=vy2vy1
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=gammaxygammav
IF (dgamma .GT. pi2) THEN
dgamma=dgammapi2
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=(dcsqs*r12*DSQRT(1dr*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=vx1m21*dvx2
vy1=vy1a*m21*dvx2
c *** velocity correction for inelastic collisions ***
vx1=(vx1vx_cm)*R + vx_cm
vy1=(vy1vy_cm)*R + vy_cm
vx2=(vx2vx_cm)*R + vx_cm
vy2=(vy2vy_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 (AZ)
CHARACTER*4 mode
INTEGER*2 error
m21=m2/m1
x21=x2x1
y21=y2y1;
vx21=vx2vx1
vy21=vy2vy1
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.0D12 should be replaced by a larger value). **********
fy21=1.0D12*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=vx1m21*dvx2
vy1=vy1a*m21*dvx2
c *** velocity correction for inelastic collisions ***
vx1=(vx1vx_cm)*R + vx_cm
vy1=(vy1vy_cm)*R + vy_cm
vx2=(vx2vx_cm)*R + vx_cm
vy2=(vy2vy_cm)*R + vy_cm
RETURN
END
