http://www.itongji.cn/cms/article/articledetails?articleid=5029

最近美国把萨德系统部署到韩国，一时心血来潮就用python模拟最简单的弹道轨迹。希望能帮助各位初学者学习python数学建模和matplotlib动态可视化模拟。

发表一下政治观点：看了战争之王的朋友可以理解，和平是军火商的噩梦。为了赚取高额军火利润，美国军火商要不停制造全球仇恨和紧张。美国在亚太不停挑拨离间各个中，日，韩，朝鲜，菲律宾，制造仇恨和冲突。这是为了能借机卖更多军火给这些国家。

暴力是人的本能之一。

洛克希德马丁公司是美国知名军火商，利润每年上百亿。萨德系统就是洛克希德马丁的产品。奥巴马和特朗普收了军火商的政治贿金，美国总统只不过是军火商的头号代理商和宣传工具。

所以不要期待民主自由的美国总统给世界带来和平，很多时候，为了赚钱，美国政客和军火商要不停制造全球冲突和仇恨。

洛克希德马丁

数学建模要用导数知识：

感谢英国大神牛顿和德国大神莱布尼茨的导数求最值方法，当导弹的瞬时速度为0时，导弹高度达到最高值（峰值），看不懂的可以去补补微积分知识，高中课本就能看懂。

Python导入math模块，表示飞行时间t_flight：

t_flight = 2*u*math.sin(theta_radians)/g

这是代码运行的界面

运行后可以观察弹道数据，设置不同发射速度和角度可以得到不同结果。

生成的动态图：

生成动态图需要导入matplotlib模块。

说明此语句意思animation.FuncAnimation(fig, update,generate,interval=5)

animation.FuncAnimation函数用于生成动态图片。fig是生成的图表对象，generate函数生成数据后传递给update函数更新，这样数据不断更新，图形也不停变化。

interval表示时间间隔，设置的值越小，运动速度越快。

代码运行平台：

Canopy python 2.7，Windows32位系统

代码汇总

源代码添加详细注解，方便各位朋友阅读理解

# -*- coding: utf-8 -*-

'''

Animate the trajectory of an object in projectile motion

'''

#seaborn增强背景效果

import seaborn

from matplotlib import pyplot as plt

from matplotlib import animation

from matplotlib.font_manager import FontProperties

import math

g = 9.8

fig = plt.figure()

ax= fig.add_subplot(111)

ax.set_aspect('equal')

#中文字体路径 设置

font=FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=14)

#获取一个列表，有205个间隔数据，每个数据间隔0.005

def get_intervals(u, theta):

intervals = []

start = 0

interval = 0.005

while start < t_flight:

intervals.append(start)

start = start + interval

return intervals

#更新时间间隔参数，从而不断改变圆的圆心坐标位置，让其移动

def update(t):

x = u*math.cos(theta_radians)*t

y = u*math.sin(theta_radians)*t - 0.5*g*t*t

circle.center = x, y

return circle,

#产生时间间隔参数，（从0,0.005,0.01一直到1.02 ）依次传递给updata函数

def generate():

for t in intervals:

yield t

def Print():

print u"初始速度（米/秒）:",u

print u"发射角度（度）",theta

print u"飞行总时间（秒）",t_flight

print u"飞行距离（米）",xmax

#初始参数，u为初始速度，theta为发射角度

u = 30

theta =60

#返回一个角度的弧度值

theta_radians = math.radians(theta)

'''

Out[65]: 0.5235987755982988

'''

#导弹飞行总时间，运用导数知识可以求得公式

t_flight = 2*u*math.sin(theta_radians)/g

intervals = get_intervals(u, theta_radians)

'''

[0,

0.005,

0.01,

0.015,

0.02,

0.025,

0.10500000000000002,

0.11000000000000003,

0.11500000000000003,

.......

0.9900000000000008,

0.9950000000000008,

1.0000000000000007,

1.0050000000000006,

1.0100000000000005,

1.0150000000000003,

1.0200000000000002]

len(intervals)

Out[67]: 205

'''

xmin = 0

#x横轴最大距离

xmax = u*math.cos(theta_radians)*intervals[-1]

ymin = 0

t_max = u*math.sin(theta_radians)/g

#y横轴最大距离

#ymax = u*math.sin(theta)*t_max - 0.5*g*t_max**2

ymax =xmax

#设置坐标轴的x,y取值范围

ax = plt.axes(xlim=(xmin, xmax), ylim=(ymin, ymax))

#创建一个圆，圆点在（0,0），半径为0.2

circle = plt.Circle((xmin, ymin), 2)

ax.add_patch(circle)

#动画函数，让炮弹不断变化，generate产生数据传递给update更新

anim = animation.FuncAnimation(fig, update,generate,interval=5)

plt.title(u'导弹发射轨迹',fontproperties=font)

plt.xlabel(u'水平距离(米)',fontproperties=font)

plt.ylabel(u'导弹运行高度（米）',fontproperties=font)

plt.show()

#输出详细参数信息

Print()

# -*- coding: utf-8 -*-
"""
Created on Wed Dec 14 09:59:34 2016
#轰炸萨德
@author: Administrator
"""
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.font_manager import FontProperties
font=FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=20)
# 创建一个fig对象，自定义fig的size
fig = plt.figure(figsize=(60,60))
# 划分fig并且选择一个子图给ax变量
ax = fig.add_subplot(1,1,1)
#width=2000,height=2000
m = Basemap(projection='mill', llcrnrlat=30, urcrnrlat=50, llcrnrlon=90, urcrnrlon=150)
m.drawcoastlines()
m.drawcountries(linewidth=2)
#m.drawrivers()
# bjlat, bjlon are lat/lon of Bei jing北京的经纬度
bjlat = 40; bjlon = 116
#tokyolat,tokyolon 表示萨德部署地的经纬度
THAADlat,THAADlon=36.119485,128.3445734
# draw parallels
m.drawparallels(np.arange(10,90,20),labels=[1,1,0,1])
# draw meridians
m.drawmeridians(np.arange(-180,180,30),labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
# fill continents, set lake color same as ocean color.
#m.fillcontinents(color='coral',lake_color='aqua')
def Draw_position(lon,lat,city,mark,markersize=100):xpt,ypt=m(lon,lat)#convert back to lat/lonlonpt,latpt=m(xpt,ypt,inverse=True)m.plot(xpt,ypt,mark,markersize) #plot a blue dot thereplt.text(xpt+100000,ypt+100000,city)
#绘制萨德坐标
Draw_position(128.3445734,36.119485,"THAAD",'c*',100)
#绘制北京坐标
Draw_position(116,40,"Beijign",'g^',100)
#链接北京和萨德的路线
m.drawgreatcircle(bjlon,bjlat,THAADlon,THAADlat,linewidth=2,color='b')
m.etopo()
#添加图例，文字说明
plt.legend(loc=4)
plt.title("轰炸萨德，制作人Toby!",fontproperties=font)
plt.show()