# The 6th homework

Chapter 2.2 problem 2.10 : cannon shell

---

---

Abstract

task01

Exercise 2.10 (L02) : Generalize the program developd for the previous problem so that it can deal with situations in which the target is at a differnt altitude than the cannon . Consider cases in which the target is higher and lower than then cannon.Also investigate how the minimum firing velocity required to hit the target varies as the altitude of the target is varied.
note add the effect of wind and give an algorithm to attack the target precisely （$θ±2°,v±5%$,$v_{wind}10%$）

task02

Develop "Super Auxiliary Precision Strike System"

---

Background

We investigate the effects of both air drag and the reduced air density at high altitudes in order to reproduce the results in textbook for L1. We generalize it to cover targets of different altitude and use Python to design the programs, which can realize the purpose of this assignment. This article deals with the ideas for these problems, the programs ,and the results.
Since the cannon shell could reach considerable high altitudes, we must be aware of the reduced air density. In this problem, we solved both the isothermal and adiabatic cases, and generalize the problem to cover targets at different altitudes.
On top of that, we can use several numerical approaches including the simple Euler method and the Runge-Kutta method to yield the numerical solutions to these problems.

---

Content

●Idea

●we assumed that the drag force for cannon shell is given by

$$F_{drag}=-B_2v^2$$
●according to Newton's second law,
$$\frac {d^2x} {dt^2}=\frac {F_{drag}cosθ} {m}=\frac {F_{drag}} {m}\frac {v_x} {v}=-\frac {B_2} {m}vv_x$$
$$\frac {d^2y} {dt^2}=\frac {F_{drag}sinθ} {m}-g=\frac {F_{drag}} {m}\frac {v_y} {v}-g=-\frac {B_2} {m}vv_y-g$$
●We also take the reduced air density at high altitudes into account. There are two mainstream models regard of this problem in statistical mechanics: isothermal ideal gas, and adiabatic approximation.
For isothermal ideal gas, the pressure p depends on altitudes according to
$$p(y)=p(0)e^{\frac {-mgy} {k_BT}}$$
●For an ideal gas the pressure is proportional to the destiny,so this leads to
$$ρ=ρ_0e^{\frac {-y} {y_0}}$$
where $y_0=\frac {k_BT} {mg}≈1.0×10^4m$ and $ρ_0$is the density at sea level≈1.225$Kg/m^3$.
●For adiabatic approximation, the density depends on altitudes according to
$$ρ=ρ_0({1-\frac {ay} {T_0}})^α$$
where a≈$6.5×10^{-3}K/m \;and\; α≈2.5$ for air. Here $T_0$ is the sea level temperature (in K)≈300K.
●Since the drag force due to air resistance is proportional to the density, so for both models we have
$$F_{drag}^{*}=\frac {ρ} {ρ_0}F_drag(y=0)$$
●Now we can get for isothermal ideal gas, the equations of motion become
$$\frac {d^2x} {dt^2}=-\frac {B_2} {m}vv_xe^{\frac {-y} {y_0}}$$
$$\frac {d^2y} {dt^2}=-g-\frac {B_2} {m}vv_ye^{\frac {-y} {y_0}}$$
●Also, for adiabatic approximation, the equations of motion become
$$\frac {d^2x} {dt^2}=-\frac {B_2} {m}vv_x({1-\frac {ay} {T_0}})^α$$
$$\frac {d^2y} {dt^2}=-g-\frac {B_2} {m}vv_y({1-\frac {ay} {T_0}})^α$$
●On top of that, we can use several numerical approaches including the simple Euler method and the Runge-Kutta method to yield the numerical solutions to these problems.

Code

# -*- coding: utf-8 -*-
"""
Created on Sun Oct 23 2016
@author: 胡胜勇
"""
import time
import math
import matplotlib.pyplot as plt
def correct(string,y_t):           ##### delete all the data of (x,y) where y<0
    x_record = string[0][-1] 
    y_record = y_t
    while True:
        if (string[1][-1] < y_t):
            if (string[1][-1] < y_t and string[1][-2] > y_t):
                x_record = string[0][-1] 
                y_record = string[1][-1]               
            del string[0][-1]
            del string[1][-1]
        else:
            string[0].append(string[0][-1]+(y_t - string[1][-1])*(string[0][-1] - x_record)/(string[1][-2] - y_record))
            string[1].append(y_t)
            break
    return string[0],string[1]
def calculate(v,theta):
    vx = []
    vy = []
    vx.append( v * math.cos(theta*math.pi/180))
    vy.append( v * math.sin(theta*math.pi/180))
    dt = 0.1
    t = []
    x = []
    y = []
    t.append(0)
    x.append(0)
    y.append(0)
    i = 1
    while (y[-1] >= 0):
        x.append(x[i - 1] + vx[i - 1]*dt)
        vx.append(vx[i - 1] - dt*0.00004*math.sqrt(vx[i - 1]**2+vy[i - 1]**2)*vx[i - 1]*(1-(0.0065*y[i - 1])/280)**2.5)
        y.append(y[i - 1] + vy[i - 1]*dt)
        vy.append(vy[i - 1]-9.8*dt -dt*0.00004*math.sqrt(vx[i - 1]**2+vy[i - 1]**2)*vy[i - 1]*(1-(0.0065*y[i - 1])/280)**2.5)
        t.append(t[i - 1] + dt)
        i+= 1
    r = y[-2] / y[-1]
    x[-1] = (x[-2] + r * x[-1]) / (r + 1)
    y[-1] = 0
    return x, y
def find_maxheight(string):
    for i in range(len(string[1])):
        if ( string[1][0] < string[1][i] ):
            string[1][0] = string[1][i]
    return string[1][0]
def scan_angle(v,theta, x_t, y_t, degree):
    ran_a = [20, 5, 2, 0.8, 0.2]
    delta = [1, 0.5, 0.1, 0.05, 0.01]
    theta = theta - ran_a[degree]
    theta_record = []
    x = []
    data = [[],[]]
    for i in range(int(ran_a[degree]*2/delta[degree]+1)):
        data = calculate(v,theta)[:]
        if ( find_maxheight(data) < y_t):
            theta = theta + delta[degree]
            x.append(100000000)
            theta_record.append(theta)
        else:
            data = correct(data,y_t)[:]
            x.append(data[0][-1])
            theta_record.append(theta)
            theta = theta + delta[degree]
    for j in range(len(x)):
        if ( abs(x[0] - x_t) > abs(x[j] - x_t)):
            x[0] = x[j]
            theta_record[0] = theta_record[j]
    return theta_record[0],x[0]   ##### debug
def scan_v(v,theta, x_target, y_target, degree):
    ran_v = [100, 10, 2, 0.8, 0.2]
    delta = [10, 2, 0.5, 0.1, 0.02]
    v = v - ran_v[degree]
    x = []
    theta_record = []
    v_record = []
    for i in range(int(ran_v[degree]*2/delta[degree]+1)):
        record_v = scan_angle(v,theta, x_target, y_target, degree)[:]
        x.append(record_v[1])
        theta_record.append(record_v[0])
        v_record.append(v)
        v = v + delta[degree]
    for j in range(len(x)):
        if ( abs(x[0] - x_target) > abs(x[j] - x_target)):
            x[0] = x[j]
            v_record[0] = v_record[j]
            theta_record[0] = theta_record[j]
    print v_record
    print theta_record
    print x
    return v_record[0],theta_record[0]
def deep_scan(v, theta, x_target, y_target, precision):
    degree = 0
    record = [[],[]]
    for i in range(precision):
        record = scan_v(v,theta,x_target,y_target,degree)[:]
        v = record[0]
        theta = record[1]
        degree = degree + 1
    return record
def judge_hitting(string,x_t,y_t):
    alpha = 0
    for i in range(len(string[0])):
        if ( x_t - 8 < string[0][-1] < x_t + 8 ):
            alpha = 1
    if (alpha == 1):
        print 'Successfully hitted'
        return 1
    else:
        print 'Miss'
        return 0
def main():
    start = time.clock()
    x_target = 14000
    y_target = 3000
#    theta0 = 180*math.atan(y_target/x_target)/math.pi
    data_record = deep_scan(700,60,x_target,y_target,5)   ### 3 for grade 1, 4 for grade 2 and 5 for grade 3
    v = data_record[0]
    theta = data_record[1]
    cannon_record = correct(calculate(v, theta),y_target)
    x_max = cannon_record[0][-1]
    judge_hitting(cannon_record,x_target,y_target)
    print '%f ' % x_max
    print '%f ' % theta
    print '%f ' %v
    end = time.clock()
    print "read: %f s" % (end - start)
    plt.figure(figsize = (8,6))
    plt.title('Trajectory of cannon shell')
    plt.xlabel('x(m)')
    plt.ylabel('y(m)')
    plt.plot(x_target,y_target,'k*',linewidth = 10,label='Target')
    plt.plot(cannon_record[0],cannon_record[1],label= 'Trajectroy')
    plt.legend(loc = 'upper left')
    plt.savefig('chapter2.png',dpi = 144)
    plt.show()
main()

Result

Target：（10000,5000）
Hitting point：（9999.244188，5000）
Angle:78.000000
Speed: 709.360000 m/s
time: 9.604472 s
figure（10000,5000）

Target：（20000,8000）
Hitting point：（20000.037368 ，8000）
Angle：47.270000
Speed:752.300000 m/s
time: 8.441148 s
figure（20000,8000）

Target：（30000,6000）
Miss
Hitting point：（ 25785.946814，6000）
Angle：50.020000
Speed:813.000000 m/s
time:read: 8.739701 s
figure（30000,6000）

precise guidance can be Accomplished in short distance(x<25000)，There is a large error when the target is in long distance.

Acknowledgments

Reference to the idea of my last homework.