Python小练习:追踪导弹仿真
警告:浏览器可能不支持html嵌入标签,那么很抱歉什么视频也看不到。建议使用最新稳定版的firefox/chromium
仿真的时候面向对象会很方便。
缘起
某天终于没有雾霾的时候,我在操场上追赶正在散步的父亲……然后,我就想多了= =跟踪导弹是个什么轨迹呢?
仿真
首先把问题简化下,如果从圆心开始追赶圆上匀速运动的物体,是什么情况。
我先自己设法用微分笔算了算,发现实在搞不定。
上网查查导弹问题看到一些简单的直线问题,都涉及一堆微分方程和欧拉法迭代啥的……
干脆自己仿真下吧。这是最初的版本,完全没有面向对象概念。
前半部分调整图像的代码完全可以不看,从while循环开始即可。
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)
我发现如果我速度够慢,未必追得上,甚至连被追踪物的轨道都不会进入……这挺出乎意料的,本来以为一定追的上
回到最初的问题,我从我的位置上追我父亲(好傻,都不知道估计下位置……)
面向对象
问题来了,如果我要仿真不只一个追踪导弹,比如还想仿真一个拦截导弹呢?
拦截失败演示
还可以用上面的方法不断扩充代码,每个对象写重复的代码。
但这时面向对象就能发挥威力,减少代码重用了。
以下是对四蜗牛聚合线问题的仿真。
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)。
搞着玩的,也许我该想想复试怎么办了……