第三章 预测物理

3.5 预测物理(Predictiong Physics)

预测子弹轨迹

一个受重力影响的子弹的轨道是一个抛物线，用公式1描述为:

$$P_t = P_0 + US_mt + \frac{1}{2}gt^2$$

其中$P_t$：时间t时子弹位置；

$P_0$：子弹初始位置；

$U$：射击是的朝向(即初始速度方向)，单位向量；

$S_m$：初始速度常量,$t$是当前时间；

$t$：当前时间；

$g$：重力加速度向量；

一般来说$g = \begin{bmatrix} 0 \ -9.81 \ 0 \end{bmatrix}\tag{1} ms^-2$，不过在游戏中跟可能会比这个值更大一些。

一元二次公式求解定理:

$$ax^2 + bx + c = 0 {\rightarrow} x = \frac{-b \pm \sqrt{b^2-4ac}}{2a}$$

根据上边的公式，可以推理出子弹到达目标高度$P_yt$的时间$t_i$:

$$t_i = \frac{-U_yS_m \pm \sqrt{U_y^2S_m^2 - 2g_y(P_{y0} - P_{yt})}}{g_y}$$

如果$t_i$没有值，即$U_y^2S_m^2 < 2g_y\sqrt{P_{0y} - P_{iy}}$,子弹无法到达位置$P_{iy}$,如图示:

如果&t_i&有一个值，即$U_y^2S_m^2 = 2g_y\sqrt{P_{0y} - P_{iy}}$,子弹只有一次经过该位置的时间,如图:

如果$t_i$有两个值，则有一次从上到下到达和一次由下到上到达,如图示：

假设目标位置为$P_y$,判断子弹是否能够到达这个位置，根据以下公式1,令

$$P_y - P_0 = \Delta$$ ,则：

$g^2t_i^4 - 4(g\Delta + S_m^2)t_i^2 + 4\Delta^2 = 0$, 即：

$$t_i = 2\sqrt{\frac{g\Delta + S_m^2 \pm \sqrt{(g\Delta + S_m^2) - g^2\Delta^2}}{2g^2}}$$

如果$g\Delta + S_m^2 < g^2\Delta^2$说明子弹不能以给于的速度够到达目标点；

如果$g\Delta + S_m^2 > g^2\Delta^2$说明有两个时间，可以到达目标点，如下图:

由于$U = \frac{2\Delta - gt_i^2}{2S_mt_i}$ 可以计算出子弹射出速度。

最终，计算子弹是否能够到达目标的实现如下代码所示:

def calculateFiringSolution(start, end, muzzle_v, gravity):
    delta = end - start

    a = gravity * gravity
    b = -4 *(gravity * delta + muzzle_v * muzzle_v)
    c = 4 * delta * delta

    # 检查是否没有结果
    if 4 * a * a > b * b:
        return None

    # 计算两个可能的时间
    time0 = sqrt((-b + sqrt(b*b - 4 * a * c)) / ( 2 * a) )
    time0 = sqrt((-b - sqrt(b*b - 4 * a * c)) / ( 2 * a) )

    # 找到相对最小的时间(最快抵达目标)
    if time0 < 0:
        if time1 < 0:
            return None
        else:
            ttt = time1
    else:
        if time1 < 0:
            ttt = time0
        else:
            ttt = min(time0, time1)

    # 返回子弹射出速度向量
    return (2*delta - gravity * ttt * ttt)/(2* muzzle_v * ttt)

弹丸阻力

增加空气阻力会使弹道的计算更加复杂。路径将不是一个抛物线，而是如下图效果：

令阻力D为：

$$D = -kv - cv^2$$

其中k是粘性系数，c是空气动力系数， v是炮弹的速度。

合并起来就是：

$${P_t}'' = g - k{P_t}' - c{P_t}'\mid {P_t}' \mid$$

第三个参数关联阻力到从一个方向到另一个方向（The third term relates drag in one direction to that in another）.

一个方向的运动不再独立于另一个方向(Motion in one direction is no longer independent of that in another)。这些导致了情况的复杂性。

有两种处理他的方式：

1. 一种解决方式是丢掉第三个变量：
   $${P_t}'' = g - k{P_t}'$$

${P_t}$ 可以被如下计算:

$${P_t} = \frac{gt - Ae^{-kt}}{k} + B$$

其中：

$$A = SmU - \frac{g}{k}$$

$B = P_0 - \frac{A}{k}$

*依据上述公式，当t=0时，$P_t$并不等于$P_0$ ，是不是有问题？应该是 $B = P_0 + \frac{A}{k}$

1. 另一种可供替换的方式是使用迭代计算，这种策略是：先使用初始的不考虑阻力的射击方程,计算出到达目标的结果，即射击角度；使用有阻力的方法根据上一步计算出的角度计算着陆点据目标点距离；重新调整这个角度；重复直到距离目标点足够近。如图示：

如果目标不在高处，那么45度是所能射出的最远距离，也是需要调整的最大角度；而最小角度是90度，此时只有朝上的速度（以我的理解是这样，但是0度是不是也是一种最小角度？虽然能射出的高度为0，这仍需要实现出来看效果）。

这种迭代计算方式的实现代码如下:

def refineTargeting(source, target, muzzleVelocity, gravity, margin):
    deltaPosition = target - source

    # 计算不考虑阻力的射击解决方案
    direction = calculateFiringSolution(source, target, muzzleVelocity, gravity)
    # 计算出此时的射击角度,由于有阻力，这个角度必定炮弹到不了目标点,所以是最小值
    minBound = asin(direction.y / direction.length())
    # 计算使用这个角度有阻力时射击能到达目标多近
    distance = distanceToTarget(direction, source, target, muzzleVelocity)
    # 检查是否已经到达目标范围
    if distance * distance < margin * margin:
        return direction
    else if distance > 0:                   # 原文是 minBoundDistance >０：
        # 炮弹超过目标，把这个角度设为最大值
        maxBound = minBound
        minBound = 0                       # 原文是minBound = -90, 我觉得如果是高射范围是在45-90，低射范围在0-45
    else:         # 炮弹未射过目标
        maxBound = 45

        direction = convertToDirection(deltaPosition, angle)    # angle是哪一个？,在这里应该是maxBound吧
        distance = distanceToTarget(direction, source, target, muzzleVelocity)

        if distance * distance < margin * margin:
            return direction
        else if distance < 0: # 无结果,可能距离太远无法抵达
            return None
    # 现在已经有了最大角度和最小角度，使用二分搜索慢慢趋向结果
    # distance = margin 这个没有用吧

    while distance * distance > margin * margin:    # 原文是distance * distance > margin * margin
        # 取角度中间值
        angle = (maxBound - minBound) * 0.5
        # 计算距离
        direction = convertToDirection(deltaPosition, angle)
        distance = distanceToTarget(direction, source, target, muzzleVelocity)

        if distance < 0:
            minBound = angle
        else:
            maxBound = angle

    return direction

def convertToDirection(deltaPosition, angle):
    # 计算平面方向(x,z)
    direction = deltaPosition
    direction.y = 0
    direction.normalize()

    # 计算垂直度数
    direction *= cos(angle)
    direction.y = sin(angle)

    return direction

依据具有阻力炮弹位置计算公式，可以计算出时间t时炮弹坐标：

def calculateDragPosition(startPos, gravity, startVelocityDir, startMuzzle_v, time):
    K = 0.15
    E = 2.718281828

    a = startPos + gravity * time / K
    b = startVelocityDir * startMuzzle_v - gravity / K
    c = (pow(E, -K * time) - 1) / K;

    return a - b * c

---

留下的问题

原书实现代码里并没有distanceToTarget的实现代码，依据阻力炮弹位置计算公式，需要先计算出炮弹高度y与目标高度y一致时的时间t1，然后计算出此时炮弹的位置以及与目标的距离，难以计算的是t1值，由于有 $power(e, -kt)$ 这种式子暂时无法根据公式计算，需要找到另外一个公式。

def distanceToTarget( start, end, direction, muzzle_v )
    # TODO 计算出$P_y = P_{targetY}$ 时的时间t1

    P = calculateDragPosition( direction, muzzle_v, t1 )

    return (P - start).length() - (end - start).length()

本节内容unity3d（版本5.6.0f3)实现demo以及c#代码下载在此，阻力炮弹预测暂未完整实现。