1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
using Plots # for plotting trajectory
# simulation parameters
Δt = 0.01 # time step
y_min = 0.0
θ_to_use = [45, 35] # in degrees
v_0 = 700.0 # in m/s
# constants
B_ref_over_m = 4.0 * 10^(-5) # in m-1, at 300K
T_ref = 300.0 # in kelvin
T_0 = 300.0 # in kelvin
g = 9.8 # in m/s^2
# isothermic parameters
y_o = 10 * 10^4 # k_BT/mg in meter
# adiabatic parameters
α = 2.5 # for air
a = 6.5 * 10^(-3) # in kelvin/meter
function adiabatic!(
x::Vector{Float64}, y::Vector{Float64},
v_y::Vector{Float64}, v_x::Vector{Float64},
t::Vector{Float64})
while y[end] >= y_min
# decompose previous positions and velocities
x_i = x[end]
y_i = y[end]
v_x_i = v_x[end]
v_y_i = v_y[end]
# calculate new positions
x_new = x_i + v_x_i * Δt
y_new = y_i + v_y_i * Δt
# calculate drag force
v_i = sqrt(v_x_i^2 + v_y_i^2)
F_drag = - B_ref_over_m * (T_0 / T_ref)^α * # temperature variation
(1 - ((a * y_i) / T_0))^α # density/altitude variation
F_drag_x = F_drag * v_x_i * v_i
F_drag_y = F_drag * v_y_i * v_i
# calculate new velocities
v_x_new = v_x_i + F_drag_x * Δt
v_y_new = v_y_i + F_drag_y * Δt - g * Δt
# store new positions and velocities
push!(x, x_new)
push!(y, y_new)
push!(v_x, v_x_new)
push!(v_y, v_y_new)
push!(t, t[end] + Δt)
end
end
function nodensity!(
x::Vector{Float64}, y::Vector{Float64},
v_y::Vector{Float64}, v_x::Vector{Float64},
t::Vector{Float64})
while y[end] >= y_min
# decompose previous positions and velocities
x_i = x[end]
y_i = y[end]
v_x_i = v_x[end]
v_y_i = v_y[end]
# calculate new positions
x_new = x_i + v_x_i * Δt
y_new = y_i + v_y_i * Δt
# calculate drag force
v_i = sqrt(v_x_i^2 + v_y_i^2)
F_drag = - B_ref_over_m # coefficient of drag alone
F_drag_x = F_drag * v_x_i * v_i
F_drag_y = F_drag * v_y_i * v_i
# calculate new velocities
v_x_new = v_x_i + F_drag_x * Δt
v_y_new = v_y_i + F_drag_y * Δt - g * Δt
# store new positions and velocities
push!(x, x_new)
push!(y, y_new)
push!(v_x, v_x_new)
push!(v_y, v_y_new)
push!(t, t[end] + Δt)
end
end
# interpolate the last point that's underground
function interpolate!(x::Vector{Float64}, y::Vector{Float64})
if y[end] == 0
return # no nothing if y is perfectly on 0
end
# calculate x_l, the interpolated x value at y=0
r = -y[end-1] / y[end]
x_l = (x[end-1] + r * x[end]) / (1 + r)
# set final values in the array to interpolated point on ground (y=0 )
x[end] = x_l
y[end] = 0.0
end
# setup empty plot to add to
p = plot(xlabel="x (m)", ylabel="y (m)", title="Cannon Shell Trajectory", xlim=(0, 30000), xticks=0:5000:30000, legend=:topright, lw=2)
for θ in θ_to_use
# arrays to store the trajectory
x = [0.0]
y = [0.0]
v_x = [v_0 * cosd(θ)]
v_y = [v_0 * sind(θ)]
t = [0.0]
# run the simulation
adiabatic!(x, y, v_y, v_x, t)
interpolate!(x, y)
plot_label = "adiabatic, θ = $θ"
plot!(x, y, label=plot_label, lw=2)
# reset arrays
x = [0.0]
y = [0.0]
v_x = [v_0 * cosd(θ)]
v_y = [v_0 * sind(θ)]
t = [0.0]
nodensity!(x, y, v_y, v_x, t)
interpolate!(x, y)
plot_label = "nodensity, θ = $θ"
plot!(x, y, label=plot_label, lw=2, linestyle=:dash)
end
# display the plot
display(p)
|
