1 files changed, 106 insertions, 18 deletions
diff --git a/balls.cl b/balls.cl
index e3450db..05a223b 100644
--- a/balls.cl
+++ b/balls.cl
@@ -1,24 +1,15 @@
 #define G_FACTOR 3600.0
 #define G (9.81 / G_FACTOR)
+#define DENSITY 1500.0
 
-float
-min(float a, float b) {
-	if (a < b)
-		return a;
-	return b;
-}
-
-float
-max(float a, float b) {
-	if (a > b)
-		return a;
-	return b;
-}
-
-float
-clamp(float x, float lo, float hi) {
-	return min(hi, max(x, lo));
-}
+float mass(float radius);
+int isCollision(float2 p1, float r1, float2 p2, float r2);
+void setPosition(float2 *p1, float r1, float2 *p2, float r2);
+float2 reaction(float2 p1, float2 v1, float m1, float2 p2, float2 v2, float m2);
+float2 unitNorm(float2 v);
+float fdot(float2 a, float2 b);
+float len(float2 v);
+float volume(float radius);
 
 __kernel void
 move(__global float2 *positions, __global float2 *velocities) {
@@ -63,6 +54,44 @@ collideWalls(__global float2 *positions, __global float2 *velocities, __global f
 }
 
 __kernel void
+collideBalls(
+	__global size_t *ballIndices,
+	__global float2 *positions,
+	__global float2 *velocities,
+	__global float *radii
+) {
+	size_t id, i1, i2;
+	float2 p1, p2, v1, v2;
+	float r1, r2, m1, m2;
+
+	id = get_global_id(0);
+	i1 = ballIndices[id];
+	i2 = ballIndices[id+1];
+
+	p1 = positions[i1];
+	p2 = positions[i2];
+	v1 = velocities[i1];
+	v2 = velocities[i2];
+	r1 = radii[i1];
+	r2 = radii[i2];
+	m1 = mass(r1);
+	m2 = mass(r2);
+
+	if (!isCollision(p1, r1, p2, r2))
+		return;
+
+	setPosition(&p1, r1, &p2, r2);
+
+	v1 = reaction(p1, v1, m1, p2, v2, m2);
+	v2 = reaction(p2, v2, m2, p1, v1, m1);
+
+	positions[i1] = p1;
+	positions[i2] = p2;
+	velocities[i1] = v1;
+	velocities[i2] = v2;
+}
+
+__kernel void
 genVertices(__global float2 *positions, __global float *radii, __global float2 *vertices) {
 	size_t ball, nsegs;
 	float2 center;
@@ -80,3 +109,62 @@ genVertices(__global float2 *positions, __global float *radii, __global float2 *
 
 	vertices[ball*get_local_size(0)] = center;
 }
+
+float
+mass(float radius) {
+	return volume(radius) * DENSITY;
+}
+
+/* Return true if the two balls are colliding. */
+int
+isCollision(float2 p1, float r1, float2 p2, float r2) {
+	float2 dist;
+	float rhs;
+
+	dist = p1 - p2;
+	rhs = r1 + r2;
+	return (dist.x*dist.x + dist.y*dist.y) <= rhs*rhs;
+}
+
+/* Set the positions of two balls at the moment of collision. */
+void
+setPosition(float2 *p1, float r1, float2 *p2, float r2) {
+	float2 mid, n;
+
+	mid = (*p1 + *p2) / 2.0f;
+	n = unitNorm(*p2 - *p1);
+	*p1 = mid - n*r1;
+	*p2 = mid + n*r2;
+}
+
+/* Return the velocity of ball 1 after colliding with ball 2. */
+float2
+reaction(float2 p1, float2 v1, float m1, float2 p2, float2 v2, float m2) {
+	float mrat, coef;
+	float2 dist;
+
+	mrat = 2.0 * m2 / (m1 + m2);
+	dist = p1 - p2;
+	coef = fdot(v1-v2, dist) / (len(dist)*len(dist));
+	return v1 - dist*mrat*coef;
+}
+
+float2
+unitNorm(float2 v) {
+	return v / len(v);
+}
+
+float
+fdot(float2 a, float2 b) {
+	return a.x*b.x + a.y*b.y;
+}
+
+float
+len(float2 v) {
+	return sqrt(v.x*v.x + v.y*v.y);
+}
+
+float
+volume(float radius) {
+	return 4.0 * M_PI_F * radius*radius*radius / 3.0;
+}