盒子
盒子
文章目录
  1. Python小练习:追踪导弹仿真
    1. 缘起
    2. 仿真
    3. 面向对象
    4. 总结

Python小练习:追踪导弹仿真

Python小练习:追踪导弹仿真

警告:浏览器可能不支持html嵌入标签,那么很抱歉什么视频也看不到。建议使用最新稳定版的firefox/chromium

仿真的时候面向对象会很方便。

缘起

某天终于没有雾霾的时候,我在操场上追赶正在散步的父亲……然后,我就想多了= =跟踪导弹是个什么轨迹呢?

操场示意图

仿真

首先把问题简化下,如果从圆心开始追赶圆上匀速运动的物体,是什么情况。

我先自己设法用微分笔算了算,发现实在搞不定。

上网查查导弹问题看到一些简单的直线问题,都涉及一堆微分方程和欧拉法迭代啥的……

干脆自己仿真下吧。这是最初的版本,完全没有面向对象概念。

前半部分调整图像的代码完全可以不看,从while循环开始即可。

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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
import matplotlib.pyplot as plt
import numpy as np
tolerance = 1e-1
radius = np.pi
v_o = 20
x_o, y_o = 0, radius

x_m, y_m = -radius, 0
v_m = 5

plt.figure(figsize=(10, 10), dpi=80)
plt.title(" missile flight simulator ", fontsize=40)
plt.xlim(-4, 4)
plt.ylim(-4, 4)
#plt.xticks([])
#plt.yticks([])

# set spines
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data', 0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data', 0))
plt.xticks([-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi], [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$+\pi/2$', r'$+\pi$'])
plt.yticks([-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi], [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$+\pi/2$', r'$+\pi$'])

# Note object and missile
plt.annotate('object start point', xy=(x_o, y_o), xycoords='data',
xytext=(+15, +15), textcoords='offset points', fontsize=12,
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.annotate('missile start point', xy=(x_m, y_m), xycoords='data',
xytext=(+15, +15), textcoords='offset points', fontsize=12,
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))

# alpha labels
for label in ax.get_xticklabels() + ax.get_yticklabels():
label.set_fontsize(16)
label.set_bbox(dict(facecolor='white', edgecolor='None', alpha=0.65))


while True:
if x_o == 0 and y_o == radius:
beta = 0
elif x_o == 0 and y_o == radius:
beta = np.pi
elif x_o < 0:
beta = np.pi / 2 * 3 - np.arctan(y_o / x_o)
else:
beta = np.pi / 2 - np.arctan(y_o / x_o)
if np.sqrt((x_o - x_m) ** 2 + (y_o - y_m) ** 2) < tolerance:
print "collision"
plt.plot(x_m, y_m, 'o')
plt.annotate('crash point', xy=(x_m, y_m), xycoords='data',
xytext=(+15, +15), textcoords='offset points', fontsize=12,
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.pause(0.1)
break
elif x_o < x_m:
alpha = np.pi + np.arctan((y_o - y_m) / (x_o - x_m))
elif x_o == x_m:
alpha = np.pi / 2
else:
alpha = np.arctan((y_o - y_m) / (x_o - x_m))
x_o = radius * np.sin(beta + v_o * 0.01 / np.pi / 2)
y_o = radius * np.cos(beta + v_o * 0.01 / np.pi / 2)
x_m = x_m + v_m * 0.01 * np.cos(alpha)
y_m = y_m + v_m * 0.01 * np.sin(alpha)
#print alpha, beta
plt.plot(x_o, y_o, 'r.', alpha=.5)
plt.plot(x_m, y_m, 'bx', alpha=.5)
plt.legend(("target", "missile"), loc="upper left", prop={'size': 12})
plt.pause(0.1)

我发现如果我速度够慢,未必追得上,甚至连被追踪物的轨道都不会进入……这挺出乎意料的,本来以为一定追的上

回到最初的问题,我从我的位置上追我父亲(好傻,都不知道估计下位置……)

面向对象

问题来了,如果我要仿真不只一个追踪导弹,比如还想仿真一个拦截导弹呢?

拦截失败演示

还可以用上面的方法不断扩充代码,每个对象写重复的代码。

但这时面向对象就能发挥威力,减少代码重用了。

以下是对四蜗牛聚合线问题的仿真。

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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
import matplotlib.pyplot as plt
import numpy as np
tolerance = 1e-1
radius = np.pi

# missile 1
x_m1, y_m1 = -np.pi, 0
v_m1 = 5

# missile 2
x_m2, y_m2 = 0, np.pi
v_m2 = v_m1
# missile 3
x_m3, y_m3 = np.pi, 0
v_m3 = v_m1
# missile 4
x_m4, y_m4 = 0, -np.pi
v_m4 = v_m1

plt.figure(figsize=(10, 10), dpi=80)
plt.title(" missile flight simulator ", fontsize=40)
plt.xlim(-4, 4)
plt.ylim(-4, 4)
#plt.xticks([])
#plt.yticks([])

# set spines
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data', 0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data', 0))
plt.xticks([-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi], [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$+\pi/2$', r'$+\pi$'])
plt.yticks([-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi], [r'$-\pi$', r'$-\pi/2$', r'$0$', r'$+\pi/2$', r'$+\pi$'])

plt.annotate('missile start point', xy=(x_m1, y_m1), xycoords='data',
xytext=(+15, +15), textcoords='offset points', fontsize=12,
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))

# alpha labels
for label in ax.get_xticklabels() + ax.get_yticklabels():
label.set_fontsize(16)
label.set_bbox(dict(facecolor='white', edgecolor='None', alpha=0.65))


class ob(object):
"""docstring for ob"""
def __init__(self, x, y):
self.x = x
self.y = y


class missile(ob):
"""docstring for missile"""
def __init__(self, x, y):
super(missile, self).__init__(x, y)

def forward(self, v, target):
"""docstring for forward"""
if self.x < target.x:
alpha = np.arctan((target.y - self.y) / (target.x - self.x))
elif self.x > target.x:
alpha = np.pi + np.arctan((target.y - self.y) / (target.x - self.x))
elif self.x == target.x and self.y < target.y:
alpha = np.pi / 2
else:
alpha = -np.pi / 2
self.x = self.x + v * 0.01 * np.cos(alpha)
self.y = self.y + v * 0.01 * np.sin(alpha)
return self.x, self.y

def distance(self, target):
"""docstring for distance"""
return np.sqrt((self.x - target.x) ** 2 + (self.y - target.y) ** 2)


class target(ob):
"""docstring for target"""
def __init__(self, x, y):
super(target, self).__init__(x, y)

def newposition(self, x, y):
"""docstring for newposition"""
self.x = x
self.y = y

m1 = missile(x_m1, y_m1)
m2 = missile(x_m2, y_m2)
m3 = missile(x_m3, y_m3)
m4 = missile(x_m4, y_m4)

while True:
if m1.distance(m2) < tolerance or m1.distance(m3) < tolerance or m1.distance(m4) < tolerance:
print "collision"
plt.plot(x_m1, y_m1, 'o')
plt.annotate('crash point', xy=(x_m1, y_m1), xycoords='data',
xytext=(+15, +15), textcoords='offset points', fontsize=12,
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.pause(0.1)
plt.show()
break
elif m3.distance(m2) < tolerance or m3.distance(m4) < tolerance:
print "collision"
plt.plot(x_m3, y_m3, 'o')
plt.annotate('crash point', xy=(x_m3, y_m3), xycoords='data',
xytext=(+15, +15), textcoords='offset points', fontsize=12,
arrowprops=dict(arrowstyle="->", connectionstyle="arc3,rad=.2"))
plt.pause(0.1)
plt.show
break
x_m1, y_m1 = m1.forward(v_m1, m2)
x_m2, y_m2 = m2.forward(v_m2, m3)
x_m3, y_m3 = m3.forward(v_m3, m4)
x_m4, y_m4 = m4.forward(v_m4, m1)
#print alpha, beta
plt.plot(x_m1, y_m1, 'bx', alpha=.5)
plt.plot(x_m2, y_m2, 'k*', alpha=.5)
plt.plot(x_m3, y_m3, 'r.', alpha=.5)
plt.plot(x_m4, y_m4, 'gp', alpha=.5)
plt.legend(("missile1", "missile2", "missile3", "missile4"), loc="upper left", prop={'size': 12})
plt.pause(0.1)

总结

面向对象方法对仿真问题非常合适,能有效简化代码,做到DRY(Don’t repeat yourself)。

搞着玩的,也许我该想想复试怎么办了……