本帖最后由 1364847132 于 2020-3-16 01:05 编辑
前言
这篇文章笔者在读初中之时就想写，奈何当时知识水平不足，只能草草放弃。如今时隔8年再来回味这个问题。
一个错误的模型
运动模型
弹弹堂的物理模型其实就是一个抛物模型。
假设没有空气阻力、不同武器发射弹丸质量相同，弹丸在飞行过程中只受重力g和风力。
由于弹丸质量相同，故风力给予弹丸的加速度与风力大小成正比。弹丸初速度与力度成正比
故
考虑最后的式子，需要拟合的参数有三个：Fn,Fw,g。
参数拟合
找到某公式计算器。可以得到许多组F,W,dx,θ。最后利用最小二乘法拟合得到Fn,Fw,g。 数据生成
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
参数拟合
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])))
跑出来的结果是:
| 参数名 |
值 |
| Fn |
0.70772385 |
| Fw |
3.92255802 |
| g |
121.764639 |
给定不同的初值p0，最后拟合出的结果不同。用不同拟合结果带入calF计算力度的结果又相同，经过验证。存在g ∝ Fn^2,g ∝ Fw。关系，所以直接给定g=10000。得到结果，
| 参数名 |
值 |
| Fn |
6.41361161 |
| Fw |
322.14261571 |
| g |
10000 |
最后可以得到力度的计算公式
带空气阻尼的运动模型
事实证明之前推断出的模型是错误的。在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个：空气阻尼系数、风力系数、重力系数。
参数拟合
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拟合来求解
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。问题就出现在这个±5°。从右往左打时，是正常的；而从左往右打时，由于会将角度进行一个转换±5°就变成了∓5°，所以造成了轨迹的不一致。下图是模拟的从右往左打和从左往右打轨迹
|
[复制链接]
发表于 2020-3-16 01:08
