22 2D中的抛物体
上一章我们在有阻力和没有阻力的条件下，建立了物体一维运动的模型。现在让我们发展到二维模型，加入棒球进行建模。
在本章中，我们考虑空气阻力建立棒球飞行模型。在接下来的章节中，我们会用这个模型去解决一系列优化问题。
22.1 棒球
为了建立棒球飞行模型，我们需要做一些建模决策。首先，我们忽略球的自旋和马格努斯力（详见 http://modsimpy.com/magnus）。在这个假设下，球在一个垂直平面内运动，也因此我们能在一个二维空间进行模拟，而不是三维空间。
空气阻力对大多数空中的抛物运动有着显著的影响，因此我们会引入一个阻力。
建立空气阻力模型，我们需要棒球质量，迎风面积和棒球的阻力系数。球的质量和直径容易得到（详见 http://modsimpy.com/baseball)。阻力系数就比较困难了，参考The Physics of Baseball（阿戴尔，The Physics of Baseball,第三版，2002年），棒球的阻力系数大概是0.33（没有单位）。
然而，这个值与速度有关。在速度较低的时候可能会高达0.5，在速度较高时可能会低至0.28。不过，这些状态的转化通常发生在我们感兴趣的速度范围内，在20m/s到40m/s之间。
尽管如此，我们从一个简单模型开始，阻力系数与速度无关。在本章末的一个练习中，你会有机会建立一个全面详尽的模型，然后看看结果与简单的模型有什么不同。
首先我们需要一个计算工具，向量对象。
22.2 向量
因为我们是在二维平面进行计算，向量会是一个很有用的计算工具。向量是一个既有大小又有方向的量。我们会用向量来表示二维或三维空间中的位置、速度、加速度和力。
ModSim库中提供一种向量对象表示向量。一个向量对象就像一个Numpy库；它包含能代表向量各个部分的元素。例如，在一个代表空间位置的向量中，属性就是x轴和y轴（如果是三维就有z轴）。一个向量对象也可能有单位，例如我们前几章举到的一些物理量。
你可以通过确定向量的组成部分来创建一个向量。下面这个向量表示一个点距右边（东）3m,距上边（北）4m：
A = Vector(3, 4) * m
你可以通过点操作来获取向量的属性：例如，A.x或A.y。例如，你也可以用括号索引的获取方式：例如，A[0]或A[1]。
类似的，你可以用点操作获取向量的大小和角度，A.mag和A.angle。向量的大小是指向量的模，如果向量表示位置，那么他的大小就是距原点的距离；如果它代表速度，那么它的大小就是速度的大小，也就是速度的值，不考虑方向。
向量的角度就是方向，表示为距x轴正方向的弧度制角。在笛卡尔平面中，弧度0就是正东方向，弧度π就是正西方向。
向量支持大多数数学运算，包括加减运算。
B = Vector(1, 2) * m
A + B
A - B
想获取这些运算的更多的定义和图解解释，可参考 http://modsimpy.com/vecops。
当你进行向量的加减运算时，ModSim库会用Numpy和Pint来检查这些向量是否具有相同的属性数和单位，这一章的笔记会展示几个向量运算的例子。
关于使用角度的一点注记：在数学中，我们用弧度来表示角度，大多数Python的功能运算也需要使用弧度。但人们通常会自然地用角度的思维方式思考，在一个项目中同时使用两种计量方式会导致不习惯，容易出错。好在Pint可以创建角度单位量来表示角度。
举个例子，从UNITS中获取一个角度，然后创建一个角度单位量来表示45度角：
degree = UNITS.degree
angle = 45 * degree
如果需要转换为弧度，可以使用to函数
radian = UNITS.radian
rads = angle.to(radian)
如果给你一个角度和速度，你可以做一个向量，用pol2cart函数从极坐标转换为笛卡尔坐标。为了证明这一点，我将分析A的角度和大小：
mag = A.mag
angle = A.angle
然后用相同的分量生成一个新的向量:
x, y = pol2cart(angle, mag)
Vector(x, y)
另一种表示A方向的方法是单位向量，它是一个向量与A指向同一方向的大小为1。你可以通过原向量除以原向量的大小来计算单位向量：
A / A.mag
我们可以用hat函数做同样的事情，因为unit矢量通常用帽子装饰，就像这样: Aˆ。
A.hat()
现在我们回到建模。
22.3 模拟棒球飞行
让我们模拟一个棒球在角度45，初始速度40 m/s。用本垒板的中心作为原点，x轴与地面平行；y轴垂直。初始高度大概是1m。
就像在21.2节中，我将创建一个Params对象，它包含系统的各项参数：
t_end = 10 * s
dt = t_end / 100
params = Params(x = 0 * m,
y = 1 * m,
g = 9.8 * m/s**2,
mass = 145e-3 * kg,
diameter = 73e-3 * m,
rho = 1.2 * kg/m**3,
C_d = 0.33,
angle = 45 * degree,
velocity = 40 * m / s,
t_end=t_end, dt=dt)
棒球的质量、直径和阻力系数都是来源第22.1节。重力加速度g是一个众所周知的量，并且空气密度rho基于海平面20℃（详见http://modsimpy.com/tempress）。我选择了t_end 的值来运行足够长的时间模拟球落地。
下面的函数使用Params对象来生成一个系统对象。这个两步过程使代码更具可读性，并使其更易于使用一些函数例如root_bisect函数。
def make_system(params):
angle, velocity = params.angle, params.velocity
# convert angle to degrees
theta = np.deg2rad(angle)
# compute x and y components of velocity
vx, vy = pol2cart(theta, velocity)
# make the initial state
R = Vector(params.x, params.y)
V = Vector(vx, vy)
init = State(R=R, V=V)
# compute area from diameter
diameter = params.diameter
area = np.pi * (diameter/2)**2
return System(params, init=init, area=area)
make_system函数用np.deg2rad将角度转换为弧度，用pol2cart计算初始速度的x和y分量。
然后它使向量对象把物体的初始位置表示为R、初始速度表示为V， 这些向量作为状态变量存储在init中。make_system还计算面积，然后用params加init和area中的所有变量创建一个系统对象。
接下来我们需要一个计算阻力的函数：
def drag_force(V, system):
rho, C_d, area = system.rho, system.C_d, system.area
mag = -rho * V.mag**2 * C_d * area / 2
direction = V.hat()
f_drag = direction * mag
return f_drag
此函数将V作为向量，并将f_drag作为向量返回。它使用阻力方程来计算阻力大小、用帽子函数来计算方向。-V.hat() 计算在V的相反方向的单位向量。
现在我们准备好斜率函数：
def slope_func(state, t, system):
R, V = state
mass, g = system.mass, system.g
a_drag = drag_force(V, system) / mass
a_grav = Vector(0, -g)
A = a_grav + a_drag
return V, A
通常，斜率函数的参数是状态对象、时间和系统对象。在这个例子中，我们不使用t，但是我们不能忽略它因为当run_ode_solver调用slope函数时，它总是提供同样的参数，不管是否需要。
State对象包含两个表示位置的状态变量R和V速度作为矢量物体。
斜率函数的返回值是这些量的导数。位置的导数是速度，所以第一个返回值是V，这是我们从状态对象中提取的。速度的导数是加速度，这就是我们要计算的。棒球的总加速度是重力和阻力引起的加速度之和。这些量既有大小又有方向，所以它们用矢量对象表示。
我们已经看到了 a_drag是如何计算的。a_grav是一个大小为重力加速度g，方向为指向负y方向的向量。使用向量来表示力和加速度使代码简洁，可读性强，不易出错。特别是当我们加上一个重力 a_grav和一个阻力a_drag时，方向很可能是正确的，因为它们是在矢量中编码的物体。而且单位肯定是正确的，否则Pint就会报告错误.
一如既往，我们可以通过在初始条件下运行来测试斜率函数：
slope_func(system.init, 0, system)
我们可以使用事件函数来停止模拟当球落地的时候。
def event_func(state, t, system):
R, V = state
return R.y
event函数使用与slope函数相同的参数，并返回R的y坐标。当y坐标通过0时模拟停止。
现在我们准备运行模拟：
results, details = run_ode_solver(system, slope_func,
events=event_func)
results是一个TimeFrame对象，每个状态变量有一列，R和V。
我们可以这样得到飞行时间：
flight_time = get_last_label(results) * s
最后的x坐标是这样的：
R_final = get_last_value(results.R)
x_dist = R_final.x
图22.1：模拟棒球飞行，位置的x和y分量时间的作用。
22.4 轨迹
我们可以这样绘制位置的x和y分量：
x = results.R.extract(✬x✬)
y = results.R.extract(✬y✬)
x.plot()
y.plot()
extract函数循环遍历R中的向量对象并从每个坐标中提取一个坐标。结果是一个时间序列。
图22.1显示了结果。正如预期的那样，x分量随着球离开本垒。y位置先上升，然后下降，在5.0 s附近下降至0 m。
查看相同数据的另一种方法是在x轴上绘制x分量以及y轴上的y分量，所以绘制的线是沿着轨迹的球穿过平面：
图22.2：模拟棒球飞行，轨迹图。
x = results.R.extract(✬x✬)
y = results.R.extract(✬y✬)
plot(x, y, label=✬trajectory✬)
图22.2显示了可视化结果的方法，即所谓的轨迹图（参见http://modsimpy.com/trajec）。
轨迹图比时间序列图更容易解释，因为它显示了炮弹的运动是什么样子的（至少从一个视角）。两个图都有用，但不要把它们弄混。如果你当你看到一个时间序列图并把它解释为一个轨迹，你就会很困惑。
在你继续之前，你可能想看看本章的笔记，第22章，ipynb，做练习。有关下载的说明运行代码，见0.4节。
本书的中文翻译由南开大学医学院智能医学工程专业2018级、2019级的师生完成，方便后续学生学习《Python仿真建模》课程。翻译人员（排名不分前后）：薛淏源、金钰、张雯、张莹睿、赵子雨、李翀、慕振墺、许靖云、李文硕、尹瀛寰、沈纪辰、迪力木拉、樊旭波、商嘉文、赵旭、连煦、杨永新、樊一诺、刘志鑫、彭子豪、马碧婷、吴晓玲、常智星、陈俊帆、高胜寒、韩志恒、刘天翔、张艺潇、刘畅。
整理校订由刘畅完成，如果您发现有翻译不当或者错误，请邮件联系changliu@nankai.edu.cn。
本书中文版不用于商业用途，供大家自由使用。
未经允许，请勿转载。
