No.6

摘要：炮弹与炮弹不在同一海拔，海拔变化自然对轨迹有影响。同时探究目标海拔改变时怎样能用最小速度打到目标，最后尝试精确击打，并计算相应误差。

背景

炮弹在空中飞行受到重力与空气阻力、风阻的影响

$$\frac{d^2x}{dt^2}=0$$ $$\frac{d^2y}{dt^2}=-g$$ $$F_{drag,x}=-B_2\mid \overline{v}-v_{wind}\mid(v_x-v_{wind})$$ $$F_{drag,y}=-B_2\mid \overline{v}-v_{wind}\mid v_y$$
通过这几个式子可解得： $x_{i+1}=x_i+v_{x,i}\Delta t$
$v_{x,i+1}=v_{x,i}-\frac{B_2\mid \overline{v}-v_{wind}\mid(v_{x,i}-v_{wind})}{m}\Delta t$
$y_{i+1}=y_i+v_{y,i}\Delta t$
$v_{y,i+1}=v_{y,i}-\frac{B_2\mid \overline{v}-v_{wind}\mid v_{y,i}}{m}\Delta t-g\Delta t$
又空气密度会随高度改变
$\rho=\rho_0(1-\frac{ay}{T_0})^\alpha$
有$F_{drag}^*=\frac{\rho}{\rho_0}F_{drag}(y=0)$
将$F_{drag}^*$以上公式，即可进行编辑代码。

正文

• 此算法中需要取定风速，空中风速随海拔变化较大，且随着地形温度等变化风速和海拔的关系式也在变化。因此此处取$v_{wind}=25m/s$,time_step=0.01, initial_v=700, T0=273，初始海拔为1000m.

• 部分代码

  import math
  import pylab as pl
  class cannon_shell1:
      def __init__(self, time_step=0.01 ,angle=1, initial_v=700, T0=273, b=0.0065, Bm=0.00004,g=9.8,c=2.5,v_wind=25):
          self.y = [1000]
          self.x = [0]
          self.t = [0]
          self.ang = angle
          self.v = [initial_v]
          self.vx =[initial_v*math.cos(self.ang)]
          self.vy =[initial_v*math.sin(self.ang)]
          self.b = b
          self.Bm = Bm
          self.T0 = T0
          self.dt = time_step
          self.g=g
          self.c=c
          self.v_wind=v_wind
      def run1(self):
          _time=0
          while (self.y[-1]>0):
              self.vx.append(self.vx[-1]-self.Bm*self.dt*math.sqrt((self.vx[-1]-self.v_wind)*(self.vx[-1]-self.v_wind)+self.vy[-1]*self.vy[-1])*(self.vx[-1]-self.v_wind)*(1-self.b*self.y[-1]/self.T0)**self.c)
              self.vy.append(self.vy[-1]-self.g*self.dt-self.Bm*self.dt*math.sqrt((self.vx[-1]-self.v_wind)*(self.vx[-1]-self.v_wind)+self.vy[-1]*self.vy[-1])*self.vy[-1]*(1-self.b*self.y[-1]/self.T0)**self.c)
              self.y.append(self.y[-1]+self.vy[-1]*self.dt)
              self.x.append(self.x[-1]+self.vx[-1]*self.dt)
              self.v.append(math.sqrt(self.vx[-1]*self.vx[-1]+self.vy[-1]*self.vy[-1]))
              self.t.append(_time)
              _time+=self.dt
  

• 绘图时选择了多个发射角度：0.5、0,6、0.7、0.8、0.9、1、1.1、1.2、1.3、1.4、1.5、。图中做了两条水平黑线，代表不同的海拔高度。

• 
• 由图中可以看出炮弹的大致打击范围。
• 加入不同速度的图像，将打击范围进行对比
• 
• 从图中可以看出不同初始速度带来的打击范围的变化。具体打击时可根据目标的位置确定最小能打至该点的速度（每个初始速度对应确定的击打范围），再通过改变射击的角度变化打击到确定目标。
• 定点击打 选定点（15000m，1000m）初始速度选择v=700m/s

• 部分代码

  class cannon_shell1:
  def __init__(self, time_step=0.01 ,angle=0, initial_v=700, T0=273, b=0.0065, Bm=0.00004,c=2.5,v_wind=25):
      self.y = [0]
      self.x = [0]
      self.t = [0]
      self.ang = angle
      self.v = [initial_v]
      self.vx =[initial_v*math.cos(self.ang)]
      self.vy =[initial_v*math.sin(self.ang)]
      self.b = b
      self.Bm = Bm
      self.T0 = T0
      self.v_wind = v_wind
      self.dt = time_step
      self.g = g
      self.c=c
  def run1(self):
      _time=0
      while (self.y[-1]<1000):
          initial_v=700
          self.ang=self.ang+0.01
          self.y = [0]
          self.x = [0]
          self.t = [0]
          self.v = [initial_v]
          self.vx =[initial_v*math.sin(self.ang)]
          self.vy =[initial_v*math.sin(self.ang)]
          while (self.x[-1]<15000):
              self.vx.append(self.vx[-1]-self.Bm*self.dt*math.sqrt((self.vx[-1]-self.v_wind)*(self.vx[-1]-self.v_wind)+self.vy[-1]*self.vy[-1])*(self.vx[-1]-self.v_wind)*(1-self.b*self.y[-1]/self.T0)**self.c)
              self.vy.append(self.vy[-1]-self.g*self.dt-self.Bm*self.dt*math.sqrt((self.vx[-1]-self.v_wind)*(self.vx[-1]-self.v_wind)+self.vy[-1]*self.vy[-1])*self.vy[-1]*(1-self.b*self.y[-1]/self.T0)**self.c)
              self.y.append(self.y[-1]+self.vy[-1]*self.dt)
              self.x.append(self.x[-1]+self.vx[-1]*self.dt)
              self.v.append(math.sqrt(self.vx[-1]*self.vx[-1]+self.vy[-1]*self.vy[-1]))
              self.t.append(_time)
              _time+=self.dt
  

• 运行时出现如下错误
  因为不知具体错误是哪所以只做到了这一步。

结论

每一个初始速度对应一个打击范围，当目标位置处于打击范围向外的边界上，能达到该目标的最小初始速度即为与边界范围相应的初始速度。