+ return 8.0 * sqrt(body->m);
+}
+
+
+/**
+ * @function Body_CheckCollision
+ * @purpose Check for collision between bodies
+ * @param A - First body
+ * @param B - Second body
+ * @param system - System that the bodies belong to
+ * @param t - pointer to variable that will store the time of closest approach
+ * @returns true if collision occurred
+ */
+bool Body_CheckCollision(Body * A, Body * B, System * system, float * t)
+{
+ assert(A != NULL && B != NULL && system != NULL && A->exists && B->exists);
+
+ if (A == B)
+ return false;
+
+ // Determine time and position of closest approach
+ float x_rel[DIMENSIONS];
+ float v_rel[DIMENSIONS];
+ float coeff_t = 0.0;
+ float c = 0.0;
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ x_rel[i] = B->x_old[i] - A->x_old[i];
+ v_rel[i] = B->v_old[i] - A->v_old[i];
+ c += v_rel[i] * x_rel[i];
+ coeff_t += v_rel[i] * v_rel[i];
+ }
+
+ *t = -c / coeff_t;
+ if (*t <= system->dt/100.0 || *t > system->dt)
+ return false;
+
+ float r2 = 0.0;
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ x_rel[i] += *t * v_rel[i];
+ r2 += pow(x_rel[i], 2);
+ }
+
+ return (r2 <= pow(Body_Radius(A), 2) + pow(Body_Radius(B), 2));
+}
+
+/**
+ * @function Body_HandleCollision
+ * @purpose Handle a collision between two Bodies
+ * @param A - The first Body
+ * @param B - The second Body
+ * @param system - The System to which the bodies belong
+ * @param t - The time of closest approach
+ */
+void Body_HandleCollision(Body * A, Body * B, System * system, float t)
+{
+ fprintf(stderr, "#Collision %p %p at t %f\n", (void*)A, (void*)B, t);
+ fprintf(stderr,"# Before Collide; A velocity is %f %f, B velocity is %f %f\n", A->v[0], A->v[1], B->v[0], B->v[1]);
+ //Calculate A's new velocity and position
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ A->v[i] = A->m * A->v_old[i] + B->m * B->v_old[i] + B->m * system->cor * (B->v_old[i] - A->v_old[i]);
+ A->v[i] = A->v[i] / (A->m + B->m);
+
+ A->x[i] = A->x_old[i] + A->v_old[i] * t; //Move to the collision point at the original velocity
+ A->x[i] += A->v[i] * (system->dt - t); //Move the rest of the step at the new velocity
+
+ }
+
+ //For cor = 0, bodies always merge; B is travelling with A
+ if (system->cor == 0.0)
+ {
+ A->m += B->m; //Increase A's mass
+ Body_Destroy(B, system); //Destroy B
+ return;
+ }
+
+ //Calculate B's new velocity and position
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ B->v[i] = B->m * B->v_old[i] + A->m * A->v_old[i] + A->m * system->cor * (A->v_old[i] - B->v_old[i]);
+ B->v[i] = B->v[i] / (B->m + A->m);
+
+ B->x[i] = B->x_old[i] + B->v_old[i] * t; //Move to the collision point at the original velocity
+ B->x[i] += B->v[i] * (system->dt - t); //Move the rest of the step at the new velocity
+
+ }
+
+ //Determine whether Body B has escaped Body A; if it hasn't, merge the two bodies
+
+
+ fprintf(stderr, "# After Collide; A velocity is %f %f, B velocity is %f %f\n", A->v[0], A->v[1], B->v[0], B->v[1]);
+
+ float v2 = 0.0;
+ float r2 = 0.0;
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ v2 += pow(B->v[i], 2);
+ r2 += pow(B->x[i] - A->x[i], 2);
+ }
+
+ fprintf(stderr, "# Collision; v is %f, escape is %f and %f\n", v2, 2.0 * system->G * A->m / sqrt(r2),2.0 * system->G * B->m / sqrt(r2));
+ if (v2 > 2.0*system->G * A->m / sqrt(r2) && v2 > 2.0*system->G * B->m / sqrt(r2))
+ return;
+
+ //Inelastically merge the two bodies
+ for (unsigned i = 0; i < DIMENSIONS; ++i)
+ {
+ A->x[i] = (A->m * A->x[i] + B->m * B->x[i]) / (A->m + B->m);
+ A->v[i] = (A->m * A->v[i] + B->m * B->v[i]) / (A->m + B->m);
+ }
+ A->m += B->m;
+ Body_Destroy(B, system);
+