- vxa = state->vx[a];
- vya = state->vy[a];
- vxb = state->vx[b];
- vyb = state->vy[b];
- nx = state->px[b] - state->px[a];
- ny = state->py[b] - state->py[a];
- if (d < -0.0001 || d > 0.0001)
- {
- cdx = nx/d;
- cdy = ny/d;
- dd = state->r[a] + state->r[b] - d;
- state->px[a] -= dd*cdx; /* just move them apart */
- state->py[a] -= dd*cdy; /* no physical rationale, sorry */
- state->px[b] += dd*cdx;
- state->py[b] += dd*cdy;
- m = sqrt (state->vx[a] * state->vx[a] +
- state->vy[a] * state->vy[a]);
- if (m < -0.0001 || m > 0.0001) /* A's velocity > 0 ? */
- {
- cosam = ((cdx * state->vx[a] + cdy * state->vy[a]) *
- state->e);
- vxa -= cdx * cosam;
- vya -= cdy * cosam; /* conserve momentum */
- vxb += cdx * cosam;
- vyb += cdy * cosam;
- }
- m = sqrt (state->vx[b] *
- state->vx[b] +
- state->vy[b] *
- state->vy[b]);
- if (m < -0.0001 || m > 0.0001)
- {
- cosam = ((cdx * state->vx[b] + cdy * state->vy[b]) *
- state->e);
- vxa += cdx * cosam;
- vya += cdy * cosam;
- vxb -= cdx * cosam;
- vyb -= cdy * cosam;
- }
+ dd = state->r[a] + state->r[b] - d;
+ /* A pair of balls that have already collided in this
+ * current frame (and therefore touching each other)
+ * should not have another collision calculated, hence
+ * the fallthru if "dd ~= 0.0".
+ */
+ if ((dd < -0.01) || (dd > 0.01))
+ {
+ cdx = (state->px[b] - state->px[a]) / d;
+ cdy = (state->py[b] - state->py[a]) / d;
+
+ /* Move each ball apart from the other by half the
+ * 'collision' distance.
+ */
+ state->px[a] -= 0.5 * dd * cdx;
+ state->py[a] -= 0.5 * dd * cdy;
+ state->px[b] += 0.5 * dd * cdx;
+ state->py[b] += 0.5 * dd * cdy;
+
+ ma = state->m[a];
+ mb = state->m[b];
+ vxa = state->vx[a];
+ vya = state->vy[a];
+ vxb = state->vx[b];
+ vyb = state->vy[b];
+
+ vela = sqrt((vxa * vxa) + (vya * vya));
+ velb = sqrt((vxb * vxb) + (vyb * vyb));
+
+ vela1 = vela * ((ma - mb) / (ma + mb)) +
+ velb * ((2 * mb) / (ma + mb));
+ velb1 = vela * ((2 * ma) / (ma + mb)) +
+ velb * ((mb - ma) / (ma + mb));
+
+ vela1 *= state->e; /* "air resistance" */
+ velb1 *= state->e;
+#if 0
+ vela1 += (frand (50) - 25) / ma; /* brownian motion */
+ velb1 += (frand (50) - 25) / mb;
+#endif
+ state->vx[a] = -cdx * vela1;
+ state->vy[a] = -cdy * vela1;
+ state->vx[b] = cdx * velb1;
+ state->vy[b] = cdy * velb1;