Python小练习：追踪导弹仿真

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

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

缘起

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

仿真

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

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

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

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

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

import matplotlib.pyplot as pltimport numpy as nptolerance = 1e-1radius = np.piv_o = 20x_o, y_o = 0, radiusx_m, y_m = -radius, 0v_m = 5plt.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 spinesax = 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 missileplt.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 labelsfor 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 pltimport numpy as nptolerance = 1e-1radius = np.pi# missile 1x_m1, y_m1 = -np.pi, 0v_m1 = 5# missile 2x_m2, y_m2 = 0, np.piv_m2 = v_m1# missile 3x_m3, y_m3 = np.pi, 0v_m3 = v_m1# missile 4x_m4, y_m4 = 0, -np.piv_m4 = v_m1plt.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 spinesax = 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 labelsfor 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 = yclass 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 = ym1 = 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)。

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