Более реалистичный 2д коллайдер / More realistic 2d collider

 В предыдущих двух постах (раз, два) я реализовал 2д-коллайдер, который применяет правило "угол падения равен углу отражения", которое визуально достаточно неплохо работает (особенно когда много частиц), но все же слишком упрощенно показывает отражения после столкновений. Например, частицы всегда имеют одну и ту же скорость и меняется только направление, а хотелось бы, чтобы частицы все же обменивались импульсами более реалистично. Поэтому я решил усложнить задачу и сделать коллайдер, который реализует более сложный закон, а именно закон абсолютно упругого столкновения и закон сохранения импульса. Все формулы и пояснения можно найти здесь. Вот ключевые моменты:

И вот эти формулы:

Собственно, это вся теория, которая нам нужна. Кто-то может заметить, что все эти иллюстрации и формулы работают для одномерного пространства, но нам же нужен двумерный случай! В википедии специально для этого указано следующее:

В случае столкновения двух тел в трёхмерном пространстве векторы импульсов тел до и после столкновения лежат в одной плоскости. Вектор скорости каждого тела может быть разложен на две компоненты: одна по общей нормали поверхности сталкивающихся тел в точке контакта, а другая параллельная поверхности столкновения. Поскольку сила удара действует только по линии столкновения, компоненты скорости, векторы которых проходят по касательной к точке столкновения, не изменяются. Скорости, направленные вдоль линии столкновения, могут быть вычислены с помощью тех же уравнений, что и столкновения в одном измерении. Окончательные скорости могут быть вычислены из двух новых компонентов скоростей и будут зависеть от точки столкновения.

То есть все, что нам в итоге надо сделать, это взять ось столкновения, которая в случае двух кругов всегда проходит через их центры и спроецировать вектор движения каждого круга на эту ось - это у нас будет изменяемая компонента, которая меняется после удара в соотв. с вышеуказанными формулами. Вторая компонтента - это проекция вектора движения каждого из кругов на тангециальную составляющую, которая вычисляется путем поворота оси столкновения на 90 градусов. Вторая компонтена у нас не меняется после удара. После того, как мы вычислили изменяемую компоненту, мы умножили этот скаляр на единичный вектор вдоль оси столкновения и сложили его с тангенциальной составляющей, получив результирующий вектор движения после удара кругов. Вот картинка, которая это поясняет:

В коде это выглядит вот так:

1. void pulseColl(const Ball& b2)
2. {
3. // direction to other ball
4. Point2f toB2 = b2.pos - pos;
5. Point2f toB2Unit = toB2.unitVector();
6. // calculate dot product
7. float dotProdB = f.dotProduct(toB2Unit);
8. // calculate speed v for both balls
9. // that is projection to hit axis
10. float v1 = dotProdB;
11. float dotProdB2 = b2.f.dotProduct(toB2Unit);
12. float v2 = dotProdB2;
13. // if move projection to hit axis
14. // for both balls are directed to move away
15. // one from other then do nothing
16. if (v1 <= 0 && v2 >= 0)
17. {
18. return;
19. }
20. // mass is equal to square of 2d ball
21. float m1 = r * r;
22. float m2 = b2.r * b2.r;
23. // speed of this ball after collision
24. float v1new = (2 * m2 * v2 + v1 * (m1 - m2)) / (m1 + m2);
25. // "to" move component
26. Point2f tmpTo = toB2Unit * v1new;
27. // "tangent" move component
28. Point2f tanUnit = Mat2x2f().rot(k_PI / 2) * toB2Unit;
29. float dotProdTan1 = f.dotProduct(tanUnit);
30. Point2f tmpTan = tanUnit * dotProdTan1;
31. // new move vector
32. Point2f newF = tmpTo + tmpTan;
33. f = newF;
34. }

Сначала мы вычислем единичный вектор вдоль оси столкновения (строки 4-5), после чего вычисляем проекцию на этот вектор нашего вектора движения (строка 7). В строке 11 мы вычисляем такую же проекцию для второго круга. В строке 16 мы проверяем, что оба обьекта двигаются не друг от друга (а если двигаются друг от друга, то ничего не делаем), Такая проверка нужна, т.к. без нее, как экспериментально выяснилось, обьекты могут "сцепляться" и кружиться в такой сцепке друг вокруг друга. Далее, в строке 24 мы производим собственно вычисления нового значения скорости вдоль оси столкновения, после чего в получаем первую компоненту (осевую или радиальную) нового значения скорости. В строках 29-30 мы получаем вторую компоненту скорости (тангенциальную) и в строке 32 вычисляем новый вектор движения путем сложения осевой и тангенциальной компонент. "Осевая" компонента потому, что вдоль оси столкновения. 

Код ядра выглядит так:

1. Ball pulseColl(Ball b1, const Ball b2)
2. {
3. // direction to other ball
4. float2 p1 = (float2)(b1.pos.x, b1.pos.y);
5. float2 p2 = (float2)(b2.pos.x, b2.pos.y);
6. float2 toB2 = p2 - p1;
7. float2 toB2Unit = getUnitVector(toB2);
8. // calculate dot product
9. float2 f1 = (float2)(b1.f.x, b1.f.y);
10. float dotProdB1 = getDotProduct(f1, toB2Unit);
11. // calculate speed v for both balls
12. // that is projection to hit axis
13. float v1 = dotProdB1;
14. float2 f2 = (float2)(b2.f.x, b2.f.y);
15. float dotProdB2 = getDotProduct(f2, toB2Unit);
16. float v2 = dotProdB2;
17. // if move projection to hit axis
18. // for both balls are directed to move away
19. // one from other then do nothing
20. if (v1 <= 0 && v2 >= 0)
21. {
22. return b1;
23. }
24. // mass is equal to square of 2d ball
25. float m1 = b1.r * b1.r;
26. float m2 = b2.r * b2.r;
27. // speed of this ball after collision
28. float v1new = (2 * m2 * v2 + v1 * (m1 - m2)) / (m1 + m2);
29. // "to" move component
30. float2 tmpTo = toB2Unit * v1new;
31. // "tangent" move component
32. float2 tanUnit = rotVect(toB2Unit, M_PI_2_F);
33. float dotProdTan1 = getDotProduct(f1, tanUnit);
34. float2 tmpTan = tanUnit * dotProdTan1;
35. // new move vector
36. float2 newF = tmpTo + tmpTan;
37. b1.f.x = newF[0];
38. b1.f.y = newF[1];
39. return b1;
40. }
41.  

Тут все более-менее то же самое. 

Исходный код находится здесь.