Python小练习:追踪导弹仿真

Published on Mar 02, 2013

Python小练习：追踪导弹仿真

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

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

缘起

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

playground.png
Figure 1: 操场示意图

仿真

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

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

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

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

前半部分调整图像的代码完全可以不看，从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)。

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

Built with in copyleft@2012-2022