一、前言
本系列是我学习物理引擎相关知识的总结。虽说最初的目标是自己实现一个完整的物理引擎，但是随着学习的深入我认识到，为了让引擎更好地模拟实际，所需要的各种数学原理和处理技巧远远超出了“涉猎”的程度。因此我决定就此止步，但是还是留下了自己的实践成果和经验总结，就在这系列文章中。
本文是“2D物理引擎基础”的第一篇，将介绍物理引擎中物体所需的物理属性
二、物理引擎中的物体
在物理引擎中，模拟的物体对象被假定为刚体。刚体是如此致密坚硬的一种物体，以至于任何碰撞都无法改变其形状。同时物理引擎还常常假定对象的密度均匀，这样对象的质心将只由其几何形状决定。通过这样的假设，物理引擎将能够通过算法优化高效地进行计算模拟：
- 刚体说明物体不会改变形状，这样当物体受到力的作用时，力的作用效果就不会造成改变物体运动之外的其他影响，如质心改变。
- 密度均匀使得质心的求解更为方便，将质心位置和偏转角度作为物体属性，就可以完全表示物体上各点的位置信息
接下来所说的物理属性也是在如上约定的基础上的。
三、运动学属性
形状的表示
物理引擎中的形状用描述物体边界的属性表示。不同形状的物体通过不同的几何参数进行表示，如矩形可用长宽表示，圆可用半径长度表示，这些较为简单。更一般地我们考虑多边形，多边形可用顶点坐标表示，但是我们希望形状是位置无关的，那么多边形的顶点坐标也不能在坐标轴上随意选取。虽然 $(0, 0), (1, 0), (1, 1), (0, 1)$ 和 $(1, 1), (2, 1), (2, 2), (1, 2)$ 表示相同形状的多边形，但我们需要选择与位置无关的那一个。我们可以选择使多边形的质心在坐标原点的顶点集合，因为我们后续需要通过质心确定运动学属性。对于均质物体来说，这也就是多边形的几何中心在原点。
我们不能人为地在设置形状时就限制顶点围成的形状中心必须在原点，这对于复杂的形状来说十分困难。我们可以根据顶点集合先确定几何中心位置，再对每一顶点减去中心坐标以将物体平移到集合中心为原点的位置。
对于多边形，计算几何中心的公式为
$$ A = \frac{1}{2} \sum_{i=0}^{N-1}(x_iy_{i+1} - x_{i+1}y_i) $$ $$ c_x = \frac{1}{6A} \sum_{i=0}^{N-1}(x_i + x_{i+1})(x_iy_{i+1} - x_{i+1}y_i) $$ $$ c_y = \frac{1}{6A} \sum_{i=0}^{N-1}(y_i + y_{i+1})(x_iy_{i+1} - x_{i+1}y_i) $$其中 $A$ 为多边形面积， $(c_x, c_y)$ 即几何中心位置。需要注意的是，在这个公式中，顶点序列为逆时针，且第一个和最后一个应是同一个坐标，以暗示这是一个闭合图形并且方便公式表示。$N$ 是不重合的顶点个数/边数。
如下是求集合中心的 python 代码
def get_geometric_center(points):
area = 0
c_x = c_y = 0
n = points.shape[0] - 1 # edge num
for i in range(n):
cross = np.cross(points[i], points[i+1])
area += cross
c_x += (points[i][0] + points[i + 1][0]) * cross
c_y += (points[i][1] + points[i + 1][1]) * cross
area = area / 2
c_x = c_x / (6 * area)
c_y = c_y / (6 * area)
return np.array([c_x, c_y])
这样得到无视位置的多边形顶点集合
mass_point = get_geometric_center(points)
points = points - mass_point
位置和运动
我们需要属性来表示物体的位置，对于二维空间而言，就需要二维的坐标表示位置；另外，我们还需要表示物体旋转的角度，在二维空间只需要一个标量即可。
对于物体的运动状态，包括速度、加速度、角速度和角加速度。前两者为二维向量，后两者为标量（虽说本质上它们是垂直于二维平面的第三维度方向的向量）。
加速度和角加速度是速度的改变方式，也就是力/力矩的作用的表现。它们随力/力矩的产生而产生，随力/力矩的消失而消失，因此并不是物体的本质属性。对物理引擎而言，我们可以在一帧里添加一定的力，由这力产生加速度与角加速度。
速度和角速度是对物体位置状态改变的度量。在物理引擎的刚体模拟中，我们特指速度为质心的速度，角速度为围绕质心旋转的角速度。
物理引擎以离散的时间间隔更新物体的位置和运动，如下所示
$$ \vec{v} = \vec{v}_0 + \vec{a}\Delta t $$ $$ \vec{p} = \vec{p}_0 +\vec{v} \Delta t $$ $$ w = w_0 + \beta \Delta t $$ $$ \theta = \theta_0 + w \Delta t $$当然，这种更新是存在误差的，我们以加速度-速度-位置为例，根据
$$ v = v_0 + \int_0^{\Delta t} a \mathrm{d} t $$ $$ p = p_0 + \int_0^{\Delta t} v(t) \mathrm{d} t $$可得
$$ p = p_0 + v_0\Delta{t} + \frac{1}{2} a \Delta{t}^2 $$可是将物理引擎的更新方式中的前两个式子可得
$$ p = p_0 + v_0\Delta{t} + a \Delta{t}^2 $$两者之间存在着误差，这是由物理引擎的离散性决定的，看似也可以直接用带有二分之一的式子来更新位置，但是这也是在假定在这段时间内加速度不变的情况下得到的，在由万有引力作用的情境下，这样的结果同样不够精确。
但是物理引擎的目的本就不是十分精确的模拟，因此这样的误差可以接受。
四、动力学属性
密度、质量和转动惯量
物体的质量由物体的形状和各点的密度决定，而物体围绕质心的转动惯量又由物体的形状和质量决定。又因为物理引擎中假设物体为刚体，物体形状确定不变，因此我们只需要确定密度就可以计算出质量和转动惯量了。
当然，对于二维，我们这里说的密度是面密度（$\sigma$）。有公式
$$ M = \int_S \sigma \mathrm{d} s $$对于物理引擎假设的均质物体
$$ M = \sigma S $$特别的，对于多边形我们有面积公式，在求几何中心时也有提及。
$$ A = \frac{1}{2} \sum_{i=0}^{N-1}(x_iy_{i+1} - x_{i+1}y_i) $$对于转动惯量，情况就有些复杂了，不同的形状有不同的计算公式。一般的形状如圆形的公式较为简单，但对于形状不定的多边形而言就比较难以确定了。但是依旧有公式，这一点似乎网上很少有文章提及
$$ I = M \frac{\sum_{n=0}^{N-1} (P_{n+1} \times P_n) ((P_n · P_n) + (P_n ·P_{n+1}) + (P_{n+1} · P_{n+1}))}{6 \sum_{n=1}^N (P_{n+1} \times P_n)} $$其中 $P_i$ 为各顶点， $\times$ 为向量叉积，$\cdot$ 为向量点积，$M$ 为质量，$N$ 依旧为不相同的顶点数/边数。
如下有代码
def get_moments_of_inertia(points, mass):
p = points
n = points.shape[0] - 1
sum1 = 0
sum2 = 0
for i in range(n):
cross = np.cross(p[i + 1], p[i])
sum1 += cross * (np.dot(p[i], p[i]) + np.dot(p[i], p[i + 1]) + np.dot(p[i + 1], p[i + 1]))
sum2 += cross
return mass * (sum1 / (6 * sum2))
这样，密度、质量和转动惯量实际只需要一个量来表示。
恢复系数与摩擦系数
我们考虑一个简单的问题，两个球体质量分别为 $m_1, m_2$，碰撞前速度为 $v_1, v_2$，碰撞后速度为 $u_1, U_2$。那么我们由动量守恒和能量守恒得
$$ m_1u_1 + m_2u_2 = m_1v_1 + m_2v_2 $$ $$ \frac{1}{2}m_1v_1^2 + \frac{1}{2}m_2v_2^2 = (\frac{1}{2}m_1v_1^2 + \frac{1}{2}m_2v_2^2) $$联立可得
$$ u_1 = \frac{(m_1 - m_2)v_1 + 2m_2v_2}{m_1 + m_2} $$ $$ u_2 = \frac{(m_2 - m_1)v_2 + 2m_1v_1}{m_1 + m_2} $$这是完全弹性碰撞的情况，我们还需要考虑存在能量损失的情况。我们引入恢复系数（coefficient of restitution） $C_r = \frac{u_1 - u_2}{v_1 - v_2}$。这代表了相对速度在碰撞后的损失率。此时有
$$ u_1 = \frac{m_1v_1 + m_2v_2 + C_r m_2 (v_2 - v_1)}{m_1 + m_2} $$ $$ u_2 = \frac{m_1v_1 + m_2v_2 + C_r m_1 (v_1 - v_2)}{m_1 + m_2} $$当 $C_r = 1$ 时，退化为完全弹性碰撞时的式子，当 $C_r = 0$ 时，碰撞后速度相同，为完全非弹性碰撞的情况。$0 < C_r < 1$ 时则为非完全弹性碰撞。
恢复系数反映了碰撞时物体弹力所产生的物体运动状态改变的效果，它与物体材质有关。当两不同材质物体发生碰撞时，我们取其恢复系数为两者中的较小值
cr = min(body_a.restitution, body_b.restitution)
在 box2d 引擎中还设置了
restitutionThreshold，只有在该值速度之上才会应用恢复系数进行反弹，当然这应该只是为了保持物体状态稳定所做的限制，并不具有实际物理意义
当两个物体接触时，在切向上会有摩擦力的作用。在物理引擎中，可以忽略动摩擦系数和静摩擦系数间的微小差别，用单一的摩擦系数表示。
摩擦系数与接触面的属性有关，因此两个不同材质的接触面可能会有所不同。为了简化程序，物理引擎可能会为每个物体设置单独的摩擦系数，并通过计算得出接触面的实际摩擦系数。如在 box2d 引擎中，就使用了将两个摩擦力相乘再开方的方式得到摩擦系数。
$$ friction = \sqrt{friction1 * friction2} $$五、总结
本文章介绍了 2d 物理引擎中的对象应该具有的基本属性，最后总结一下这些属性。
- 形状：物体的几何形状，不同形状表示方式不同，其中多边形用顶点序列表示
- 位置：（质点的）空间位置和（相对质点的）旋转角度
- 速度：（质点的）线速度和（相对质点的）角速度
- 加速度：线加速度和角加速度，由力所产生
- 密度：物体的面密度
- 质量：质量由密度和形状求出
- 转动惯量：转动惯量由质量和形状求出
- 恢复系数：表示物体的弹性程度
- 摩擦系数：表示物体的粗糙程度
由此我们可以定义物理引擎的主体类 Body。为了方便，这里只使用多边形作为形状。
class Body:
force = np.zeros(2)
torque = 0
def __init__(self, points, position, density,
rotate=0, v=0, w=0, restitution=0, friction=0.66):
mass_point = get_geometric_center(points)
points = points - mass_point
self.points = points
self.position = position
self.density = density
self.mass = density * get_area(self.points)
self.inv_mass = 1 / self.mass
self.inertia = get_moments_of_inertia(self.points, self.mass)
self.inv_inertia = 1 / self.inertia
self.rotate = rotate
self.v = v
self.w = w
self.restitution = restitution
self.friction = friction
def apply_force(self, force, pos):
self.force += force
self.torque += np.cross(pos - self.position, force)
def step(self, dt):
self.v = self.v + self.force / self.mass * dt
self.position = self.position + self.v * dt
self.w = self.w + self.torque / self.inertia * dt
self.rotate = self.rotate + self.w * dt
self.force = np.zeros(2)
self.torque = 0
注意这里并未存储线加速度和角加速度，而是存储了力和力矩。