网站 › 『编程语言区』 › 浅析弹弹堂物理模型

1364847132 发表于 2020-3-16 01:08

浅析弹弹堂物理模型

本帖最后由 1364847132 于 2020-3-16 01:05 编辑

# 前言
这篇文章笔者在读初中之时就想写，奈何当时知识水平不足，只能草草放弃。如今时隔8年再来回味这个问题。
# 一个错误的模型
## 运动模型
弹弹堂的物理模型其实就是一个抛物模型。
假设没有空气阻力、不同武器发射弹丸质量相同，弹丸在飞行过程中只受重力g和风力。
由于弹丸质量相同，故风力给予弹丸的加速度与风力大小成正比。弹丸初速度与力度成正比
故

考虑最后的式子，需要拟合的参数有三个：Fn,Fw,g。
## 参数拟合
找到[某公式计算器](http://news.4399.com/baiduapp/dandantang/gongsi/)。可以得到许多组F,W,dx,θ。最后利用最小二乘法拟合得到Fn,Fw,g。
数据生成
```python
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 - f, 50


def calfix45(dx, f):
    temp =
    return temp - 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 - f, 30


def calflo65(dx, f):
    temp =
    return temp, 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)
            W.append(cf)
            A.append(temp)

    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, )# 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, f30, f45, f65, hi))
    F = np.concatenate((f50, f30, f45, f65, hi))
    W = np.concatenate((f50, f30, f45, f65, hi))
    A = np.concatenate((f50, f30, f45, f65, hi))
    return dx, F, W, A

```
参数拟合
```python
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 =
p = leastsq(error, p0, args=(dx, F, W, A))
print(p)
print(calF(p, np.array(), np.array([-3.0]), np.array()))
```
跑出来的结果是:

参数名|值
:-:|:-:
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个：空气阻尼系数、风力系数、重力系数。
## 参数拟合
```python
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, )
      assert time != 0
      return time

    times = np.array()
    return computePosition(vx, W * wind, R, times)


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


p0 =
p = leastsq(error, p0, args=(dx, F, W, A))
print(p)
```
最后得到拟合结果

参数名|值
:-:|:-:
R|0.89927083
W|5.8709153
g|-172.06527992
由于运动方程是个超越方程，所有最后力度的求解也只能通过fsolve拟合来求解
```python
    def _calF(angel, wind, dx, dy):
      R, W, g =
      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, )
                assert time != 0
                return time

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

      f = fsolve(solve, )
      if f > 100:
            raise ValueError
      return f
```
# 附言

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

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

Nineteenhundred 发表于 2020-4-18 02:05

楼主你的65度力度表是错的。30和50的非常对。65的力度表：
12        21        27        33        37        42        46        50        54        57
60        63        66        69        72        76        78        81        84        87

工地搬砖本科生 发表于 2020-4-1 23:06

1364847132 发表于 2020-3-30 10:56
我已经明确指出，这里是单方向的运功方程了。v0x，v0y，有意义吗？x方向和y方向都是这个运动方程，我为什 ...

仔细算了一遍，又和研究抛体的论文对比了下，我起初看你的公式有那么多的觉得奇怪的地方，其实只是你的角标引起误会，也怪我愚钝。
含风情况下的参数方程消参应该可以用泰勒公式把exp(x)那玩意展开几项就能把有风情况的轨迹近似解得到。（前几天出去了，来不及回复不好意思）泰勒展开有时间我得再整整

riwfhiu 发表于 2020-3-16 02:13

爷的青春回来了{:17_1085:}

saga64 发表于 2020-3-16 03:16

相当不错的文章，质量很高

aixiaodemj 发表于 2020-3-16 07:37

感谢分享~哈啊 好怀念这个游戏啊，是不是那些写瞄准器的就是靠楼主这篇文章写的

eyey000 发表于 2020-3-16 08:28

Zeaf 发表于 2020-3-16 09:28

青春啊....
啊，数学知识有限，物理就更加了。。
这种研究值得鼓励{:1_893:}

zeddit 发表于 2020-3-16 10:17

好文，分析的很细

yxdyxd163 发表于 2020-3-16 12:36

传说中理工男。。厉害

lvchao0329 发表于 2020-3-18 21:51

哈哈，不如把核弹的概率调100%{:1_927:}{:1_927:}

步步高打火机 发表于 2020-3-19 19:39

楼主可以的，厉害

查看完整版本: 浅析弹弹堂物理模型