## Fortran Program Code for 2D Elastic and Inelastic Collision of 2 Balls
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 The equations on which the code is based have been derived at c http://www.plasmaphysics.org.uk/collision2d.htm 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 |