Computational Physics 2.10L1

孔德锋（2014301020081）

一、摘要

利用欧拉公式，考虑海拔高度对迎面风阻和重力加速度的影响，并引入迎面风阻的影响，计算发射角度为40度时炮弹命中直线距离为100km，垂直距离为1km，0km，-1km的不同高度的三处所需要的最小发射速度，从而得出最小发射速度随目标的海拔高度变化而变化的规律。

二、背景

2.10：Generalize the program developed for thr previous problem so that it can deal with situations in which the target is at a different altitude than the cannon.consider cases in which the target is higher and lower than cannon.Also investigate how the minimum firing velocity required to hit the target varies as the altitude of the target is varied.

三、正文

设置初值，发射角度为40度，迎面风阻vw=-100m/s，T取值300K（地面温度）

class projectile:
    def __init__(self,h,time_step=0.1,intial_x=0,vw=100,\
    intial_y=0,intial_T=300):
        self.x=[intial_x]
        self.y=[intial_y]
        self.dt = time_step
        self.T=intial_T
        self.h=h
        self.angle=40 #发射角度
        self.vw=vw #引入迎面风阻

我们需要求出40度时打到目标大发射速度，已知目标位置，求初速度，可利用扫描的方法，扫描范围可依据无阻力时打到目标位置的初速度往上扫描为标准，因为相同发射角度，打到同一目标，有阻力时要比无阻力时需要的初速度更大，这里我们选择900m/s~1200m/s（扫描上限不宜太高，因为(1-a*self.y[-1]/self.T)**2.5会出现negative number cannot be raised to a fractional power的错误（分数幂内不能为负数)。

def v_angle(self):
        self.er=[]
        for v in range(900,1200,1): #速度扫描范围
            self.angle
            self.v_x=[v*cos(self.angle*pi/180)-self.vw]
            self.v_y=[v*sin(self.angle*pi/180)]
            while(1): #以下为引可变入阻力和变化的重力加速度的欧拉方程
                self.v=((self.v_x[-1]-self.vw)**2+self.v_y[-1]**2)**0.5
                self.x.append(self.v_x[-1]*self.dt+self.x[-1])
                self.y.append(self.v_y[-1]*self.dt+self.y[-1])
                self.g=g0*(R/(R+self.y[-1]))**2
                self.den=(1-a*self.y[-1]/self.T)**2.5
                self.F_drag_x=-B2*self.den*self.v*self.v_x[-1]
                self.F_drag_y=-B2*self.den*self.v*self.v_y[-1]
                self.v_x.append(self.v_x[-1]+self.F_drag_x*self.dt)
                self.v_y.append(self.v_y[-1]+(-self.g+self.F_drag_y)*self.dt)
                if self.x[-1]>100000:
                    self.er.append(abs(self.y[-1]-self.h))
                    break
        print 900+self.er.index(min(self.er))

目标点我们选择三处为（100000，1000），（100000，0），（100000，-1000）。

b= projectile(1000)
c= projectile(0)
d= projectile(-1000)
b.v_angle()              
c.v_angle()       
d.v_angle()

四、结果

打到目标点（100000，1000）初速度为1004m/s
（100000，0）初速度为958m/s
（100000，-1000）初速度为931m/s
从上可以看出目标高度月底所需初速度越小，这点是符合常识的，可以想象，这和炮弹打在同一高度，水平距离越小，所需初速度越小是等价的。

进一步改进

1.从结果和代码中可以看出初速度只精确到各位（其实也不精确，因为是在时长dt=0.1下得到的），因为我们有两点可改进：对for v in range(900,1200,1)可改进速度步长为0.1；当炮弹到达目标点附近时可进行精确迭代即设置时间间隔dt=0.01.
2.上述只是在给出发射角度下，求最小速度，更一般的是，我们可以不限制角度求最小发射速度。方法基本类似，只是多了一次扫描，一次对发射速度，一次对发射角度。
代码如下：
方法同上，只做两点说明：
发射角度扫描范围的选取:由于我们已经对2.6，2.7，2.8，2.9进行研究，发现在有无租赁时最远反射距离的发射角度基本在40都到60度之间，因此为运算更快，选择40`60度。
发射速度：还是以无阻力时发射速度往上扫描，因发射角度范围增大，发射速度也要相应扩大 选择（800~1600m/s）m/s

from math import cos, sin, pi
a=6.5e-3
g0=9.8
B2=4e-5
R=6.3e6
class projectile:
    def __init__(self,h,time_step=0.1,intial_x=0,\
    intial_y=0,intial_T=300):
        self.x=[intial_x]
        self.y=[intial_y]
        self.dt = time_step
        self.T=intial_T
        self.h=h
    def v_angle(self):
        self.er=[]
        for v in range(800,1600,1):#第二次扫描 对发射速度，分度为1m/s
            self.angle
            self.v_x=[v*cos(self.angle*pi/1800)]
            self.v_y=[v*sin(self.angle*pi/1800)]
            while(1):
                self.v=(self.v_x[-1]**2+self.v_y[-1]**2)**0.5
                self.x.append(self.v_x[-1]*self.dt+self.x[-1])
                self.y.append(self.v_y[-1]*self.dt+self.y[-1])
                self.g=g0*(R/(R+self.y[-1]))**2
                if 1-a*self.y[-1]/self.T>0:
                    self.den=(1-a*self.y[-1]/self.T)**2.5
                else:
                    self.den=0
                self.F_drag_x=-B2*self.den*self.v*self.v_x[-1]
                self.F_drag_y=-B2*self.den*self.v*self.v_y[-1]
                self.v_x.append(self.v_x[-1]+self.F_drag_x*self.dt)
                self.v_y.append(self.v_y[-1]+(-self.g+self.F_drag_y)*self.dt)
                if self.x[-1]>50000:
                    self.er.append(abs(self.y[-1]-self.h))
                    break
        self.v_ang=700+self.er.index(min(self.er))
    def vmin_angle(self):
        self.v=[]
        for self.angle in range(400,600,1):#第一次扫描 对发射角度 分度为0.1度
            self.v_angle()
            self.v.append(self.v_ang)
        print min(self.v)
b=projectile(1000)
c=projectile(0)
d=projectile(-1000)
b.vmin_angle()              
c.vmin_angle()       
d.vmin_angle()                    

六、感谢

老师的PPT