第六次作业

摘要

本文讨论了考虑空气阻力以及空气阻力随海拔变化的情况下（恒温模型），进行定点射击的参数设置。一门加农炮发射时的全部参数只有初始速度（包括速度的方向，对于二维情况就是速度的横纵坐标），而目标可以用海拔跟距离两个值表示。本文解决了以下问题：只要给定目标位置，即可给出全部能击中目标的速度的组合（用一条二维曲线表示出）。

代码

def accelerate(x,y,vx,vy):
    v=sqrt(vx*vx+vy*vy)
    ansax=-0.00004*v*vx
    ansay=-g-0.00004*v*vy
    return (ansax,ansay)
import pylab
from math import sqrt
from math import cos
from math import sin
from math import pi
from math import atan
g=9.8
class Euler2Dsimulation:
    def __init__(self,speed,theta,aimx,aimy,acc=accelerate,timestep=0.001):
        self._timestep=timestep
        self._x=[0]
        self._y=[0]
        self._vx=[speed*cos(theta)]
        self._vy=[speed*sin(theta)]
        self._aimx=aimx
        self._aimy=aimy
        self._accelerate=acc
    def calculate(self):#速度过小返回false
        while True:
            self._x.append(self._x[-1]+self._vx[-1]*self._timestep)
            self._y.append(self._y[-1]+self._vy[-1]*self._timestep)
            (ax,ay)=self._accelerate(self._x[-1],self._y[-1],self._vx[-1],self._vy[-1])
            self._vx.append(self._vx[-1]+self._timestep*ax)
            self._vy.append(self._vy[-1]+self._timestep*ay)
            if self._x[-1]>=self._aimx:
                return self._y[-1]-self._aimy
            if self._vy[-1]<=0 and self._y[-1]<=self._aimy:
                return -1000
    def plot(self):
        pylab.plot(self._x, self._y)
        pylab.xlabel('x ($m$)')
        pylab.ylabel('y ($m$)')
def dichotomy(para1,para2,func,eps=1):
    temp=(para1+para2)/2
    ft=func(temp)
    while abs(ft)>eps:
        f1=func(para1)
        f2=func(para2)
        if f1*ft>0 :
            para1=temp
        else:
            para2=temp
        temp=(para1+para2)/2
        ft=func(temp)
    return temp
max_speed=5000
target=(15000,2000)
min_theta=atan(target[1]/target[0])
ltheta=[]
lspeed=[]
theta=min_theta
while theta<pi/2:
    theta+=0.01
    if Euler2Dsimulation(max_speed,theta,target[0],target[1]).calculate()<0 :
        continue
    ltheta.append(theta)
    def func(para):
        return Euler2Dsimulation(para,theta,target[0],target[1]).calculate()
    lspeed.append(dichotomy(10,max_speed,func))
lx=[]
ly=[]
for i in range(len(ltheta)):
    lx.append(lspeed[i]*cos(ltheta[i]))
    ly.append(lspeed[i]*sin(ltheta[i]))
pylab.plot(lx,ly)
pylab.show()

结果展示

在这里，我们假定目标的坐标为（15000,2000），然后给出了所有可能打到这个位置的初始发射速度。图像曲线上的每一个点都代表了一种发射速度，从而可以轻易看出我们想要的信息，比如最小发射速度等。
（这里不得不吐槽python的运行速度，我跑了三十分钟都没出结果。真！的！慢！所以不得已我将精度调的非常低，让运行时间可以接受）

（这幅图上点的横纵坐标代表发射速度的横纵分量，按这些速度发射即可射中目标。理应是光滑曲线，但限于精度只取了几个离散点）