So the problem was DotProduct. It returned x1y2 + x2y1, which is nonsense. It should be x1x2 + y1*y2.
The updated code is here :
https://github.com/EdgarReynaldo/Interceptor
https://github.com/EdgarReynaldo/Interceptor/blob/master/Circle.cpp
The collision detection is super easy. Here is the code :
#SelectExpand
1double CalculateCollision
(Circle
* c1 , Circle
* c2
) {
2 if (!c1
|| !c2
) {return -1.
0;}
3
4 const double dx
= c2->cx
- c1->cx
;
5 const double dy
= c2->cy
- c1->cy
;
6 const double dvx
= c2->vx
- c1->vx
;
7 const double dvy
= c2->vy
- c1->vy
;
8 const double rsq
= (c1->rad
+ c2->rad
)*(c1->rad
+ c2->rad
);
9
10 /// Quadratic equation Ax^2 + Bx + C = 0
11 const double A
= (dvx
*dvx
+ dvy
*dvy
);
12 const double B
= 2.
0*(dvx
*dx
+ dvy
*dy
);
13 const double C
= dx
*dx
+ dy
*dy
- rsq
;
14
15 if (A
== 0.
0) {
16 return -1.
0;/// No relative movement, collision is impossible
17 }
18
19 /// Overlap check here - if they are already colliding at time 0.0, let them pass undisturbed
20 if (C
< 0.
0) {return -1.
0;}
21
22 const double TWOA
= 2.
0*A
;
23 const double DISCRIM
= B
*B
- 4.
0*A
*C
;
24 if (DISCRIM
< 0.
0) {return -1.
0;}/// no real roots
25 else if (DISCRIM
== 0.
0) {
26 return -B
/TWOA
;
27 }
28 const double SQRTD
= sqrt(DISCRIM
);
29 const double T1
= (-B
- SQRTD
)/TWOA
;
30 const double T2
= (-B
+ SQRTD
)/TWOA
;
31 if (T1
>= 0.
0) {
32 return T1
;
33 }
34/// else if (T2 >= 0.0) {
35/// return T2;
36/// }
37 return -1.
0;/// No collision possible in the future
38}
And the new code for bouncing is here :
#SelectExpand
1bool MakeCirclesBounce2(Circle* c1 , Circle* c2) {
2 if (!c1 || !c2) {return false;}
3
4 /// Normalized normal vectors, these point from one circle towards the other
5 Vec2 N1 = Vec2(c2->cx - c1->cx , c2->cy - c1->cy).Normalize();
6 Vec2 N2 = Vec2(c1->cx - c2->cx , c1->cy - c2->cy).Normalize();
7
8 Vec2 V1(c1->vx , c1->vy);
9 Vec2 V2(c2->vx , c2->vy);
10
11 const double m1 = c1->mass;
12 const double m2 = c2->mass;
13
14 Vec2 I1 = V1*m1;
15 Vec2 I2 = V2*m2;
16/// double Itotal = (I1 + I2).Magnitude();
17
18 /// The angle between V and N determines how much energy is transferred
19
20 const double I1M = I1.Magnitude();
21 const double I2M = I2.Magnitude();
22
23 Vec2 I1N(0,0);/// Momentum of circle one in the normal direction
24 Vec2 I2N(0,0);/// Momentum of circle two in the normal direction
25
26 if (I1M > 0.0) {
27 I1N = ScalarProjection(I1 , N1);/// Magnitude of N is always 1, they're normalized
28 }
29 if (I2M > 0.0) {/// Object has momentum
30 I2N = ScalarProjection(I2 , N2);/// Magnitude of N is always 1, they're normalized
31 }
32
33 if (DotProduct(I1,N1) > 0.0) {
34 /// Circle one is moving towards circle two
35 /// Give circle ones normal momentum to circle two
36 I1 = I1 - I1N;/// This momentum is lost, transferred to the other circle
37 I1N *= ELASTICITY;/// Energy is lost due to inelasticity
38 if (c2->fixed) {/// We hit an immovable object
39 I1N *= -1.0;/// Reflection of normal energy
40 I1 += I1N;/// rebound effect
41 }
42 else {
43 I2 = I2 + I1N;/// The remaining momentum is gained by the other circle
44 }
45 }
46
47 if (DotProduct(I2,N2) > 0.0) {
48 /// Circle two is moving towards circle one
49 /// Give circle twos normal momentum to circle one
50 I2 = I2 - I2N;/// This momentum is lost
51 I2N *= ELASTICITY;/// Energy lost due to inelasticity
52 if (c1->fixed) {/// We hit an immovable object
53 I2N *= -1.0;/// Reflection of normal energy
54 I2 += I2N;/// rebound effect
55 }
56 else {
57 I1 = I1 + I2N;/// The remaining momentum is gained by the other circle
58 }
59 }
60/// double Itotal2 = (I1 + I2).Magnitude();
61
62 bool changed = false;
63 if (!c1->fixed && m1 > 0.0) {
64 V1 = I1*(1.0/m1);
65 c1->SetSpeed(V1.x , V1.y);
66 changed = true;
67 }
68 if (!c2->fixed && m2 > 0.0) {
69 V2 = I2*(1.0/m2);
70 c2->SetSpeed(V2.x , V2.y);
71 changed = true;
72 }
73 return changed;
74}
The code for ScalarProjection is here :
Vec2 ScalarProjection(Vec2 A , Vec2 B) {
B.Normalize();
return B*DotProduct(A,B);
}
And the code for DotProduct is here :
inline double DotProduct(const Vec2& v1 , const Vec2& v2) {
return v1.x*v2.x + v1.y*v2.y;
}
What's really interesting though, is the collision table. Check it out :
https://github.com/EdgarReynaldo/Interceptor/blob/master/CollTable.cpp
https://github.com/EdgarReynaldo/Interceptor/blob/master/CollTable.hpp
It is a vector of pairs of circles, representing each possible combination of N circles. It is exactly (N*(N-1))/2 in size. If you're doing a 1000 circle collision resolution, it will make 499,500 pairs. But the beauty is you only have to recalculate it if the velocity of one of the circles changes. This eliminates 99% of the wasted calculations of doing a frame by frame overlap check.
EDIT
{"name":"611747","src":"\/\/djungxnpq2nug.cloudfront.net\/image\/cache\/6\/3\/635673b9af24f3c2c81797b18d2333ee.png","w":1027,"h":800,"tn":"\/\/djungxnpq2nug.cloudfront.net\/image\/cache\/6\/3\/635673b9af24f3c2c81797b18d2333ee"}
You'll notice that balls 13 and 14 are missing. They escaped. 
{"name":"611748","src":"\/\/djungxnpq2nug.cloudfront.net\/image\/cache\/c\/9\/c9d7ea71f2e32e14391e295f5548be94.png","w":1027,"h":800,"tn":"\/\/djungxnpq2nug.cloudfront.net\/image\/cache\/c\/9\/c9d7ea71f2e32e14391e295f5548be94"}
The fuller it is the more often they escape.
It's really quite stable, except for the escaping bit.
Aha! I fixed it. No more escape for you my pretties. 
{"name":"611749","src":"\/\/djungxnpq2nug.cloudfront.net\/image\/cache\/a\/a\/aab5277509ad27b99a5096f0cd412ec9.png","w":1026,"h":801,"tn":"\/\/djungxnpq2nug.cloudfront.net\/image\/cache\/a\/a\/aab5277509ad27b99a5096f0cd412ec9"}
The problem was in CollTable.cpp on line 114 :
#SelectExpand
102std::vector<CollInfo*> CollTable::GetFirstCollisionsEarlierThanDT(double dt) {
103 std::vector<CollInfo*> clist;
104 const unsigned int N = ctable.size();
105 double first = -1.0;
106 for (unsigned int n = 0 ; n < N ; ++n) {
107 CollInfo& info = ctable[n];
108 if (info.dt < 0.0) {continue;}
109 if (info.dt > dt) {continue;}
110 if (first < 0.0) {
111 first = info.dt;
112 }
113 if (info.dt > first) {
114- break;
115+ continue;
116 }
117 clist.push_back(&ctable[n]);
118 }
119 std::sort(clist.begin() , clist.end() , CompareCInfo);
120 return clist;
121}
https://github.com/EdgarReynaldo/Interceptor/commit/2e72f6c26576b87a82db4c0acbe8d992a658a051#diff-c7e697629af47a878a6163ae0d115670L114
]]>
