浅析弹弹堂物理模型

前言

这篇文章笔者在读初中之时就想写,奈何当时知识水平不足,只能草草放弃。如今时隔8年再来回味这个问题。

一个错误的模型

运动模型

弹弹堂的物理模型其实就是一个抛物模型。
假设没有空气阻力、不同武器发射弹丸质量相同,弹丸在飞行过程中只受重力g和风力。
由于弹丸质量相同,故风力给予弹丸的加速度与风力大小成正比。弹丸初速度与力度成正比

考虑最后的式子,需要拟合的参数有三个:$F_n,F_w,g$

参数拟合

找到某公式计算器。可以得到许多组$F,W,d_x,\theta$。最后利用最小二乘法拟合得到$F_n,F_w,g$。

数据生成

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
import numpy as np


def calfix50(dx, f):
temp = [0, 14.1, 20.1, 24.8, 28, 32.5, 35.9, 39.0, 42.0, 44.9, 48.3, 50.5, 53, 55.5, 58, 60.5, 63.0, 65.5, 68.0,
70.0, 72.5]
return temp[dx] - f, 50


def calfix45(dx, f):
temp = [0, 0, 0, 20, 28.7, 30, 0, 35, 0, 0, 45, 0, 0, 0, 55, 57, 0, 60, 0, 0, 68]
return temp[dx] - f, 45


def calfix30(dx, f):
temp = [0, 14, 20, 24.7, 28.7, 32.3, 35.7, 38.8, 41.8, 44.7, 47.5, 50.2, 52.8, 55.3, 57.9, 60.3, 62.7, 65.7, 67.5,
69.8, 72.1]
return temp[dx] - f, 30


def calflo65(dx, f):
temp = [0, 13, 21, 26, 31.5, 37, 41, 44, 48.5, 53, 56, 58, 61, 64, 67, 70, 73, 76, 79, 82, 85]
return temp[dx], 65 + 2 * f


def high(dx, f):
return 90, 90 - dx + 2 * f


def generate(func, distence):
dx = [] # 屏距
F = [] # 力度
W = [] # 风力
A = [] # 角度
for f in range(-59, 60):
cf = f / 10
for x in distence:
temp = func(x, cf)
dx.append(x)
F.append(temp[0])
W.append(cf)
A.append(temp[1])

return np.array(dx), np.array(F), np.array(W), np.array(A)


def gen(): # 生成数据
f50 = generate(calfix50, range(1, 21)) # 50定角
f45 = generate(calfix45, [3, 4, 5, 7, 10, 14, 15, 17, 20]) # 45定角
f30 = generate(calfix30, range(1, 21)) # 30定角
f65 = generate(calflo65, range(1, 21)) # 65变角
hi = generate(high, range(1, 21)) # 高抛

dx = np.concatenate((f50[0], f30[0], f45[0], f65[0], hi[0]))
F = np.concatenate((f50[1], f30[1], f45[1], f65[1], hi[1]))
W = np.concatenate((f50[2], f30[2], f45[2], f65[2], hi[2]))
A = np.concatenate((f50[3], f30[3], f45[3], f65[3], hi[3]))
return dx, F, W, A

参数拟合

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
from scipy.optimize import leastsq
import numpy as np
from dataGenetator import gen
import matplotlib.pyplot as plt

dx, F, W, A = gen()


def func(p, F, W, A):
Fn, Fw, g = p # 力度系数、风力系数
angel = A * np.pi / 180
sin = np.sin(angel)
cos = np.cos(angel)

V0 = Fn * F
t = 2 * V0 * sin / g
a = Fw * W
return V0 * cos * t + 0.5 * a * t * t


def calF(p, dx, W, A): # 力度计算
Fn, Fw, g = p # 力度系数、风力系数
angel = A * np.pi / 180
sin = np.sin(angel)
cos = np.cos(angel)
a = Fw * W
temp = dx / (2 * (sin * sin * a + sin * cos * g))
return g / Fn * np.sqrt(temp)


def error(p, dx, F, W, A):
return dx - func(p, F, W, A)


p0 = [10, 10, 10]
p = leastsq(error, p0, args=(dx, F, W, A))
print(p)
print(calF(p[0], np.array([16]), np.array([-3.0]), np.array([59])))

跑出来的结果是:

参数名
$F_n$ 0.70772385
$F_w$ 3.92255802
$g$ 121.764639

给定不同的初值p0,最后拟合出的结果不同。用不同拟合结果带入calF计算力度的结果又相同,经过验证。存在$g \propto F_n^2,g \propto F_w$。关系,所以直接给定$g=10000$。得到结果,

参数名
$F_n$ 6.41361161
$F_w$ 322.14261571
$g$ 10000

最后可以得到力度的计算公式

带空气阻尼的运动模型

20-3-3更新,事实证明之前推断出的模型是错误的。在60°左右的半抛时,力度过大,在80°左右的高抛时,力度过小。

运动模型

好在网上有弹弹堂泄露的服务器源码。服务器是.net写的,在包Game.Logic.Phy.Maths.EulerVector.cs中有如下函数:

m是物体质量,af是空气阻尼系数,f是物体所受外力(x方向是风力,y方向是重力),dt是时间片段,x0是物体所在位置,x1是物体速度,x2是物体加速度。
假设物体质量m=1,可以得到如下公式:

由此推导得出了物体在单方向(x、y方向)上的运动方程。物体在x方向上受持久力风力,在y方向上受持久力重力,故需要拟合的参数同样是3个:空气阻尼系数、风力系数、重力系数。

参数拟合

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
from scipy.optimize import leastsq
from scipy.optimize import fsolve
import math
import numpy as np
from dataGenetator import gen

dx, F, W, A = gen()


def func(p, F, wind, angel):
R, W, g = p # 空气阻力系数,风力系数,重力系数
angel = angel * math.pi / 180
vx = np.cos(angel) * F
vy = np.sin(angel) * F

def computePosition(v0, f, R, t):
temp = f - R * v0
ert = np.power(math.e, -R * t)
right = temp * ert + f * R * t - temp
return right / (R * R)

def getTime(v0):
solve = lambda t: computePosition(v0, g, R, t)
time = fsolve(solve, [100000])
assert time[0] != 0
return time[0]

times = np.array([getTime(i) for i in vy])
return computePosition(vx, W * wind, R, times)


def error(p, dx, F, W, A):
return dx - func(p, F, W, A)


p0 = [10, 10, -10]
p = leastsq(error, p0, args=(dx, F, W, A))
print(p)

最后得到拟合结果

参数名
$R$ 0.89927083
$W$ 5.8709153
$g$ -172.06527992

由于运动方程是个超越方程,所有最后力度的求解也只能通过fsolve拟合来求解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def _calF(angel, wind, dx, dy):
R, W, g = [0.89927083, 5.8709153, -172.06527992]
angel = angel * math.pi / 180

def solve(F):
vx = math.cos(angel) * F
vy = math.sin(angel) * F

def computePosition(v0, f, R, t):
temp = f - R * v0
ert = np.power(math.e, -R * t)
right = temp * ert + f * R * t - temp
return right / (R * R)

def getTime(v0):
solve_l = lambda t: computePosition(v0, g, R, t) - dy
time = fsolve(solve_l, [100000])
assert time[0] != 0
return time[0]

t = getTime(vy)
return computePosition(vx, W * wind, R, t) - dx

f = fsolve(solve, [100])
if f[0] > 100:
raise ValueError
return f[0]

在阅读源码的过程中,同样解释了我多年来的疑惑。为什么从左到右和从右到左的三叉轨迹不相同?

三叉弹的另外两发分别是+5°、力度乘以1.1和-5°、力度乘以0.9。问题就出现在这个$\pm5°$。从右往左打时,是正常的;而从左往右打时,由于会将角度进行一个转换$\pm5°$就变成了$\mp5°$,所以造成了轨迹的不一致。下图是模拟的从右往左打和从左往右打轨迹

文章目录
  1. 1. 前言
  2. 2. 一个错误的模型
    1. 2.1. 运动模型
    2. 2.2. 参数拟合
  3. 3. 带空气阻尼的运动模型
    1. 3.1. 运动模型
    2. 3.2. 参数拟合